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1.  INTRODUCTION 


A.  OVERVIEW 

The  analysis  of  electromagnetic  problems  in  the  time  domain  has  become  more 
common  in  recent  decades,  hi  this  approach,  the  differential  or  integral  equations  for  the 
particular  problem  under  consideration  are  typically  solved  numericalfy  u^g  a  short 
time- duration  waveform  as  excitation.  By  far  the  most  common  of  these  techniques  is  the 
so-called  finite  difference  time-domam  (FD-TD)  method. 

In  the  FD-TD  method,  the  problem  is  spatially  discretized  ("sanqiled”)  using  a  spatial 
step  Al,  usually  in  a  Cartesian  coordinate  system  and  tenqiorally  “sarrqiled"  u^g  a  time 
step  At.  The  differential  operators  in  the  differential  form  of  Maxwell's  equations  are 
approximated  by  finite  differences  or,  equivalently,  integrals  in  the  integral  form  of 
Maxwell's  equations  are  approximated  by  finite  sums.  In  this  research,  a  pulse  waveform  is 
excited  at  the  center  of  the  discretization  grid,  and  the  fields  at  the  grid  nodes  are 
computed  at  discrete  time  steps.  The  field  at  a  particular  grid  node  at  time  t  can  be 
determined  firom  the  fields  of  this  and  adjacent  nodes  calculated  for  file  previous  time  step 
t  -  At.  The  electric  and  magnetic  fields  obtained  via  the  FD-TD  method  are  the  result  of  a 
"marching  in  time"  process  A^ilOse  accuracy  is  affected  by  the  approximations  introduced 
at  each  time  step.  These  include  the  grid  size,  time  step,  length  of  time  observation, 
excitation  waveform,  and  so  forth. 
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There  are  several  advantages  to  the  FD-TD  method.  First,  it  is  relatively  simple  to 
implement  for  complex  objects  because  of  spatial  discretization.  Inhomogeneous  media  is 
handled  by  assigning  varying  electric  and  magnetic  properties  (s,  p,  and  a)  to  individual 
cells.  Although  the  FD-TD  electric  and  magnetic  field  "update"  equations  appear  complex, 
they  contain  only  addition,  substraction,  and  multiplication.  Ftirthermore,  conq)uter 
memory  requirements  are  usually  less  demanding  than  those  for  the  method  of  mommits 
(MOM). 

An  impnrfa-nt  consideration  in  the  FD-TD  method  is  how  to  terminate  the  spatial 
discretization  grid.  The  FD-TD  equation  derived  in  the  following  section  will  describe 
how  waves  propagate  along  an  infinite  grid.  However,  in  all  practical  applications  the  grid 
will  have  to  be  restricted  to  a  finite  number  of  nodes.  Therefore,  a  method  for  grid 
"termination"  is  required.  The  "standard"  FD-TD  update  equations  are  valid  for  all  nodes 
except  those  on  the  grid  edges,  because  the  grid  edge  nodes  have  fewer  neighbors  than  the 
non-edge  (internal)  nodes.  Therefore,  edge  nodes  with  a  different  number  of  neighbors 
have  update  equations  different  fi:om  the  equations  derived  for  the  internal  grid  nodes.  The 
question  is  how  can  one  determine  the  equations  for  the  edge  nodes?  One  approach  is  to 
try  to  absorb  a  wave  incident  firom  inside  the  grid  onto  the  boundary  such  that  there  is  no 
reflection  back  into  the  grid.  This  is  referred  to  as  an  absorbing  boundary  condition 
(ABC)  [Ref  2].  However,  the  grid  itself  is  not  ideal,  that  is  a  wave  does  not  pass  through 
the  grid  without  some  reflection  even  for  those  nodes  that  are  not  at  grid's  edges  (this  is 
the  consequence  of  ^atial  and  temporal  field  sampling).  Requiring  that  there  is  no 
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reflection  off  the  grid  edges  is  thus  too  restrictive,  because  there  is  some  reflection  at  each 
and  every  node  within  the  grid.  We  thus  propose  a  different  approach  that  we  refer  to  as 
the  Tranq)arent  Grid  Termination  (TGT).  We  start  by  making  only  one  demand  on  the 
grid  termination  condition:  that,  ideally,  the  solution  within  the  grid  should  be  the  same 
solution  that  one  would  obtain  with  an  infinitely  large  grid,  or  in  other  words,  that 
reflections  off  a  grid  edge  node  are  the  same  as  the  reflections  off  any  other  node  inside 
the  grid.  Note  that  if  the  grid  were  infinitely  large,  the  grid  edges  would  be  at  infinity  and 
the  equations  for  the  grid  edge  nodes  would  be  irrelevant,  since  the  excited  or  scattered 
field  would  take  an  infinitely  long  time  to  reach  the  grid  edges,  and  thus  the  edge 
equations  would  never  really  be  needed.  The  TGT  requirement  is  conceptually  very  anq}le 
and  fairly  straightforward  to  implement.  We  will  use  the  multiport  concept  (used  in 
analysis  of  linear  systems)  to  introduce  the  TGT  concept,  and  details  vsdll  be  shown  in  the 
following  sections. 

B.  OBJECTIVE 

The  performance  of  the  TGT  in  1-D,  2-D,  and  3-D  can  be  assessed  based  on  the 
residual  reflection  off  the  TGT  boundary.  In  this  thesis  the  total  power  of  the 
electromagnetic  field,  excited  by  a  unit  an:q)Utude  delta  pulse  in  the  grid  center,  will  be 
observed  as  a  function  of  time.  If  the  grid  were  ideal,  the  energy  of  the  field  at  any  node 
will  be  zero  for  all  time  after  the  pulse  has  passed  that  node.  The  total  energy  within  an 
infinite  ideal  grid  would  be  constant  however.  The  total  energy  within  a  finite  ideal  grid 
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with  ideal  grid  termination  condition  will  be  constant  until  the  pulse  has  reached  the  grid 
edge  and  will  then  quickly  go  to  zero  as  the  pulse  leaves  the  grid.  However,  because  the 
grid  is  not  ideal,  there  are  some  residual  values  at  any  grid  node,  even  after  the  initial  pulse 
has  passed  that  node.  Therefore,  the  energy  within  a  non-ideal  grid  will  not  be  constant 
but  will  fluctuate  about  its  value  for  the  ideal  grid.  These  fluctuations  generally  decrease 
with  time  but  they  would  exist  even  if  the  grid  were  infinite.  If  the  grid  is  terminated,  the 
residual  energy  within  the  grid  will  increase  because  of  additional  reflections  off  the  grid 
edges.  The  global  effect  of  the  grid  termination  can  thus  be  assessed  by  corcparing  the 
total  grid  powers  as  fimctions  of  time  for  the  following  two  cases:  an  "infinite"  grid,  and  a 
finite  size  grid  with  TGT  or  ABC.  Note  that  the  "infinite"  grid  is  simulated  by  using  a 
large  grid  and  finishing  the  observations  before  the  pulse  that  started  at  the  grid's  cmiter 
has  reached  any  of  the  grid  edge  nodes.  Note  that  if  the  TGT  or  ABC  were  ideal,  there 
would  be  no  increase  in  energy  within  the  grid  after  the  pulse  has  reached  the  grid  edges. 
If  the  TGT  is  not  ideal,  the  relative  change  of  the  aiergy  within  the  grid  (compared  to  the 
case  of  an  "infinite"  grid)  is  a  measure  of  TGT  "quality":  the  smaller  the  change,  the  better 
the  TGT  would  be.  We  will  thus  compare  the  total  energies  within  the  TGT  boundary  for 
a  large  grid  with  Dirichlet  boundary  condition  (any  boundary  condition  is  appropriate  for 
this  "infinite"  grid  simulation  because  the  reflections  off  the  edges  of  the  large  grid  do  not 
reach  the  smaller  grid  boimdary  within  the  observation  time-window)  and  the  smaller  grid 
with  the  TGT  boundary.  Once  again,  the  large  grid  must  provide  sufficient  distance 
between  the  TGT  and  Dirichlet  boundaries  such  that  the  reflections  of  the  Dirichlet 
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boundary  do  not  reach  the  TGT  boundary  within  the  observation  interval.  The  measure  of 
residual  reflection  off  the  TGT  boundary  will  be  the  difference  of  the  energies  within  the 
TGT  boundary  for  the  "infinite"  grid  and  the  TGT.  It  is  convenient  to  normalize  this 
difference  to  the  energy  within  the  TGT  boundary  for  the  "infinite"  grid,  and  express  the 
resulting  ratio  in  dB.  The  objective  of  this  thesis  is  thus  to  determine  the  quality 
(measured  in  dB  as  described  above)  of  the  Transparent  Grid  Termination  (TGT) . 
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n.  ANALYSIS  OF  TGT  FOR  1-D  FD-TD 


A.  FD-TD  FORMULATION  IN  1-D 

1.  FD-TD  Equations  in  1-D 

The  incident  and  the  scattered  electromagaetic  fields  and  the  media  parameters  in 

1-D  problems  can  vary  with  one  spatial  coordinate  only.  We  will  denote  this  coordmate 

as  z.  The  fields  will  also  be  fimctions  of  time  t.  The  media  will  be  assumed  stationary, 

that  is,  the  media  parameters  are  not  fimctions  of  time  t.  1-D  electromagnetic  fields  must 

— >  — > 

be  Transverse  Electro-Magnetic  (TEM)  fields,  that  is  the  field  coroqponents  E  and 
must  be  in  a  plane  transverse  to  the  direction  of  field  propagation  [Ref  5].  Since  we  have 

denoted  the  direction  of  propagation  as  z,  the  fields  can  be  deaoted  as  transverse-to-z  or 

TEM^.  The  electric  and  magnetic  vectors  of  a  TEM  field  must  be  perpendicular  to  each 

other.  The  E  and  H  field  vectors  and  the  direction  of  propagation  form  a  tr^let  of 

orthogonal  vectors,  as  shown  below. 
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A  ^ 

E 


Figure  2.1.  TEM  Field  Components. 

Tlie  axes  were  selected  such  that  the  electric  field  vector  is  x-directed  and  the 
magnetic  field  vector  is  y-directed. 

E(z,t)=:E,(z,t)x  (2.1) 

(2.2) 

The  electric  and  magnetic  fields  satisfy  Maxwell's  curl  equations  vriuch  express  fiie 
electric  and  magnetic  field  coupling  and  can  be  written  in  the  integral  or  the  differential 
form  We  will  use  the  integral  forms  because  they  are  applicable  to  finite  (line,  surface,  or 
volume)  domains  wdiereas  the  differential  forms  are  applicable  to  infinitesimally  small 
domains  (points).  Shown  below  are  the  integral  forms  of  Maxwell's  curl  equations  for  a 
TEM^field. 

_ ^  r  y 

Exiz,  fyx  •  dh  =  [xiz)Hy(z,  t)^  •  dsE  ■  (2-3) 

Hy{z, t)y  •  din  =  z(z)Ex(z, t)x  •  dsn  ■  +  a(z)Ex{z, ijx  •  dsn  (2.4) 

The  second  equation  does  not  have  the  source  current  term  since  it  has  been 
assumed  that  there  were  no  source  currents  in  the  domain  of  interest.  The  contours  of 
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integration  for  the  electric  and  the  magnetic  field  contour  integrals  are,  in  general,  different 
and  are  thus  labeled  Q  for  an  electric  field  circulation  and  Cj,  for  a  magnetic  field 
circulation.  (A  closed  path  line  integral  will  be  referred  to  as  the  circulation.)  Similarly,  the 
surfaces  associated  with  the  contours  are  labeled  Sg  for  the  magnetic  field  surfece  integral 
and  Sji  for  the  electric  field  surface  integral  Surface  integrals  on  the  right-hand  side  will  be 
referred  to  as  fluxes.  The  first  step  in  discretizing  the  curl  equations  is  to  select  the 
contours  for  the  electric  and  magnetic  field  circulation.  Since  the  essence  of  Maxwell's 
curl  equations  is  the  coupling  of  the  electric  and  magnetic  fields,  the  contours  will  be 
selected  such  that  this  coupling  is  achieved  in  a  straightforward  manner,  as  shown  below. 


Figure  2.2.  Contours  for  Electric  and  Magnetic  Field  Circulation  and  Fluxes. 

The  electric  field  contoius  Cg  and  the  magnetic  field  contours  Cjj  are  square  (Alby 
Al)  contovus,  but  in  orthogonal  planes.  The  Cg  contours  are  in  the  xz-plane,  while  the  C„ 
contours  are  in  the  yz-plane.  The  contours  can  be  conopared  to  links  of  a  chain,  evoking 
the  idea  of  electric  and  magnetic  field  linkage.  This  selection  of  such  contours  leads  to  the 
so-called  dual  discretization  grid,  also  referred  to  as  the  Yee  lattice  [Ref  1]. 
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2.  Derivation  of  Magnetic  Field  Update  Equation 


The  discrete  equivalent  of  the  electric  field  drculation  will  be  determined  next.  A 
contour  is  shown  below. 


X 

A 


y  “  - 

E(z-Al/2,t)  ^  E(z-+Ai/2,t) 


Figure  2.3.  A  Contour  for  Electric  Field  Circulation  and  Magnetic  Flux  Calcxilations. 

The  center  of  the  surfece  Sg  is  assumed  to  have  the  coordinates  (x,  z).  A  local 
coordinate  system  (^,Q  can  be  establidied,  with  the  origin  at  the  center  of  the  contour. 
Any  point  (x',z')  within  or  on  the  contour  Cg  (or,  equivalently,  any  point  on  the  surfece 

Sg)  can  be  specified  by  its  local  coordinates  ^  and  C 

x'=x  +  ^  z'  =  z+C  and 

The  local  coordinates  will  be  used  in  evaluation  of  the  line  and  surface  integrals 

that  constitute  the  integral  forms  of  Maxwell's  equations.  The  electric  field  is  not  a 

function  of  the  x-coordinate  and  the  electric  field  is  zero  along  the  top  and  the  bottom 

sides  of  the  contour  because  there  can  be  no  z-diiected  field  conq)onents  of  a  TEM^  field. 

The  circixlation  of  the  electric  field  around  the  electric  field  contour  Cg  ,  assuming 

coxmter-clockwise  reference  direction  such  that  the  normal  to  the  surface  Sg  is  in  the 

— ^ 

direction  of  the  magnetic  field  //(z,  t),  can  be  evaluated  exactly 
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I  +C,  0^  •  ^  =  n  +  f  ’  ~  f  ’  0^  •  (2.5) 

•'Cif  •'-J  2 

5,(z+C,0l?  =  [^.(z+f,0-^.(^-f,0]  •  A/  (2.6) 

The  rate  of  change  of  magnetic  flux  through  the  surface  (within  the  contour  C^) 

can,  on  the  other  hand,  be  evaluated  only  approximately,  since  the  exact  way  that  the 
magnetic  flux  density 

B{z,  t)  =  Vi(z)H{z,  i)  (2.7) 

varies  with  the  z-coordinate  is  not  known  a  priori  Note  that  there  is  no  variation  of  the 

magnetic  flux  density  with  the  x-coordinate.  We  will  next  evaluate  the  magnetic  flux 

approximately 

n(z+Off,(2+C()7**r  =  r!  (2«) 

The  surface  integral  reduces  to  a  line  integral,  since  the  integrand  is  not  a  function 
of  the  local  x-coordinate  ^ 

H(z  +  0^>'(^  +  co7 •  3^  =  A/ 0  p(z  +  O^(y(^  +  C0^C  (2.9) 

There  are  infinitely  many  ways  to  model  the  variation  of  the  magnetic  flux  density 

Avith  the  local  z-coordinate  ^  within  the  contour.  Since  the  magnetic  flux  density  B(z+^,t) 
is  the  product  of  the  permeability  |j,(z+^)  and  die  magnetic  field  H(z+(^,t)  we  first  need  to 
assume  a  certain  variation  of  the  magnetic  field  with  ^  such  that  the  integral  can  be 
evaluated  over  S^.  The  sinoplest  model  assumes  that  the  contour  width  A1  is  small  enough 
such  that  the  magnetic  field  H(z+^,t)  within  flie  contour  may  be  assumed  constant  and 
equal  to  the  value  at  the  contour's  center  H(z,t).  This  is  equivalent  to  u^g  a  piece-wise 
constant  ("pulse"  expansion)  approximation  of  die  actual  magnetic  field  variation  with  the 
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z-coordinate.  The  above  assumption  allows  that  the  magnetic  field,  although  constant 
within  a  contour,  can  change  from  one  contour  to  an  other.  This  yields  an  approximate 
expression  for  the  for  the  magnetic  flux 

IT  (210) 

Sb  2 

The  integral  of  the  permeability  p(z+C)  can  be  re-written  in  the  following  manner 

J_*|^(^  +  Qd^  =  A^[iJ.t^(z+0<c]  (2.11) 

The  term  in  the  brackets  is  recognized  as  the  average  permeability  within  the  contour 
Ce- 

H«v*(r)  =  ^J;|n(z+CK  (2.12) 

The  approximate  expression  for  the  magnetic  flux  through  can  now  be  written  using  the 
average  permeabihty  as 

[f  \x.(z  +  Q)Hy(z + Cj  0  T  *  dsE  ^  i^cn>g(z')Hy(z,  t)  (2. 13) 

JJ  Sb 

Again  this  approximate  expression  resulted  from  the  piece-wise  constant 
approximation  of  the  magnetic  field  with  respect  to  the  z-coordinate.  The  main  advantage 
of  the  piece-wise  constant  ejqransion  employed  above  is  its  sin^licity.  Better  accuracy  can 
be  achieved  by  using  more  involved  models  for  the  field  variation  with  z,  but  at  the 
expense  of  increasing  the  conq)lexity  and  conqrutational  time.  The  first  curl  equation  of 
Maxwell  can  now  be  replaced  by  an  approximate  equation 

which  can  be  gimplified  (because  of  media  stationarity)  to 

Ex(z  +  y,  /)  —  Ex(z  —  y,  /)  »  ^Hy(z,  t)}  •  A/  •  (lovg(z)  (2. 15) 

This  equation  can  also  be  re-wntten  as 
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(2.16) 


The  time  derivative  operator  in  the  above  equation  is  typically  replaced  by  the 
finite  difference  approximation  [Ref  1].  However,  an  alternate  approach  is  presented  here, 
such  that  the  approximation  of  the  field  tenq)oral  variation  is  shown  to  be  analogous  to  the 
approximations  already  introduced  for  the  field  spatial  variation.  The  above  equation  is 
integrated  with  respect  to  time  to  get  an  approximate  equation  for  die  magnetic  field  at  the 
present  time  t 

+  2’  ~  ^ 

A  similar  integral  equation  can  be  written  for  the  magnetic  field  at  the  same  spatial 
location  but  at  an  earlier  time  t-At 

-  r  -  T-  (218) 

Subtracting  the  two  equations  we  get 

Hyiz,  t)  -  Hy(z,  r  -  AO  «  +  f’  “  f’  (2- 19) 

which  can  be  also  written  as  an  "update"  equation  (assuming  that  the  previous  field  value 
at  the  same  location  is  known) 

Hy(z,  t)  ^  Hy(z,  E,(z  +  f .  t)*  -  E.(z  -  f ,  t)*]  (2.20) 

The  integrals  on  the  right-hand  side  can  not  be  evaluated  exactly,  because  the 
exact  temporal  variation  of  the  electric  fields  within  the  At  interval  prior  to  t  is  generally 
not  known,  just  like  the  spatial  integral  for  the  magnetic  flux  could  not  have  been 
evaluated  exactly  because  the  exact  variation  of  the  magnetic  field  within  the  contour  was 
not  known.  However,  the  integrals  can  be  evaluated  approximately  by  assuming  a  certain 
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variation  of  the  electric  field  with  the  ten:q)oral  variable  x.  The  sim5)lest  approach, 
consistent  with  the  assurcqptions  made  for  the  field  spatial  variation,  would  be  to  assume 
that  the  interval  At  is  small  enough  such  that  the  electric  field  can  be  assumed 
approximately  constant  within  At  and  equal  to  the  value  at  the  center  of  the  time  interval 
(t-At,  t) 


£(z  +  f,T)s 

=  £(z  +  f,r-f) 

for 

t-At<x<t 

(2.21) 

£(z-f.T)> 

.£(z-f,/-f) 

for 

t-At<x<t 

(2.22) 

The  above  represents  a  piece-wise  constant  approxiniation  of  the  electric  field  variation 
with  respect  to  the  tercporal  variable  t.  The  approximate  expression  for  the  magnetic  field 


at  location  z  and  the  time  t  now  becomes 


Hy(z,  I)  »  W/z,  /  -  AC)  -  +  f ,  r  -  f )  -  £.(z  -  f .  /  -  f )]  (2.23) 

The  above  equation  can  be  written  in  a  somewhat  different  form,  using  the  identity 
1  ,  1  1  1  11  1 _ = 


—  “^TT  =  (2.24) 

,nSn  lUn  Zn  M-r  M-r 


where  the  free  space  velocity  of  propagation  is  denoted  v^,  the  intrinsic  impedance  of  free 
space  is  denoted  Zq  (Zo  «  377f2),  and  p,,  denotes  relative  permitivity.  Introducing  the  grid 
"propagation"  velocity 


_  A/ 


(2.25) 


we  can  write  the  approximate  "update"  equation  for  the  magnetic  field  as 

Hy(z,  /)  H,(z,  ( -  A/)  -  +  f ,  ( -  f )  -  £.(z  -  f ,  <  -  f )]  (2.26) 

The  equation  sinplifies  for  the  case  of  non-magnetic  media  (p =1) 

Hy{z,  t)  «  Hy(z,  t-At)-  Yo^[e,(z  +  ^,t-f)-E,(z-f,t- f )]  (2.27) 
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The  last  two  equations  show  that  the  magnetic  field  H  and  the  electric  field  E  are 
evaluated  at  points  shMed  spatially  by  Al/2,  and  at  instants  separated  in  time  by  At/2.  The 
relationship  between  q)atial  and  temporal  "san^les"  of  the  electric  and  the  magnetic  fields 
is  thus  the  same  the  samples  are  shifted  with  respect  to  each  other  by  one-half  of  the 
sanq)ling  interval  This  is  the  essence  of  the  Yee  lattice  and  appHes  to  "san:q)ling"  of  2-D 
and  3-D  fields  as  well  Dual  approximate  equations  for  the  electric  field  "updates"  can 
now  be  derived  using  the  same  procedure  shown  for  the  magnetic  field  update  equations. 

3.  Derivation  of  Electric  Field  Update  Equation 

We  start  with  the  Maxwell's  curl  equation  for  the  circulation  of  the  magnetic  field 
Hy(z, iyx  =  t)x  •dsH  +  o(z)E^(z, t)x  •  dsn  (2.28) 

The  discrete  equivalent  of  the  magnetic  field  circulation  will  be  determined  first, 
using  the  contour  shown  below. 

Ch 

/  . - . >, 

I  if  H(zHAi/2,t) 


Figure  2.4.  A  Contour  for  Magnetic  Field  Circulation  and  Electric  Fliix  Calculations. 

A  local  coordinate  ^stem  (v|/,  Q  can  be  established,  with  the  origin  at  the  carter  of 
the  contour.  Any  point  (y’,z')  within  or  on  the  contour  Cjj(or,  equivalent^,  any  point  on 

the  surfece  S^)  can  be  q)eci£led  by  its  local  coordinates 

y' =y  +  \^  z' where  -^<vj/<+^  and 


The  circulation  of  the  magnetic  field  around  ,  assuming  counter-clockwise 

reference  direction  such  that  the  normal  to  the  surface  is  in  the  direction  of  the  electric 

— ^ 

field  E  (z,  f),  is  therefore  given  exactly  by  the  following  single  expression 

i  Hy(z  +  U)y  •■d^v'y  (2.29) 

•'”2  2 

I  Hy{z  +  i:,J)y  =  (2.30) 

»  C// 

The  rate  of  change  of  the  electric  flux  through  the  surface  Sy  can  only  be  evaluated 

approximately,  since  the  exact  way  that  the  electric  flux  density 

D(z,t)  =  z{z)Hiz,t)  (2.31) 

varies  in  the  z-direction  is  not  known  a  priori.  First  we  evaluate  the  flux  integrals 

z(z+QE,(z+U)x  . 2//  =  f| s(z  +  Q£,(2  +  C,0'x  . (2.32) 

Jj  Sfj  2  2 

fj  <j(z  +  OE;,(z + Q,  fix  •dsH=Vli  cj(z + QEx{z + C,  0^  •  ~xdy}dL,  (2.33) 

Since  the  integrands  are  not  fimctions  of  the  local  y-coordinate  \|/  the  surface  integrals 
sinqilify  to  line  integrals 

IT  s(z  -I-  QExiz + tyx  •  dsfi  =  A/  f  ^  8(z  +  OEx(z + C,  f)dl^  (2.34) 

JJSh  •'“2 

IT  a(z  -I-  OEx(z + C,  0^  •dsn  =  Al  cr(z + QExiz  +  C,  t)dt;  (2.35) 

JJ  Sff  ^  -J 

Next,  we  assume  that  the  contour  width  A1  is  small  enough  such  that  the  electric 
field  E(z+C,t)  within  the  contour  may  be  assumed  constant  and  equal  to  the  value  at  the 
contour's  center  E(z,t).  This  assumption  yields  approximate  expressions  for  the  for  the 
flux  integrals 

ff  8(z  +  QEx{z + (yx  •  dsH  =  A/  -  Ex{z,  t)  •  \_l  8(z + QdC,  (2.36) 
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(2.37) 


o(z + 0£.(z + C. =  A/ •  0  ■  |l|  <!(z + CK 

or,  using  the  average  permitivity  and  conductivity 

^^^z{z+QEx{z+C,,t)x  •  dsH  =  (Atf  ■Zavg{z)-Ex(z,t)  (2.38) 

C5(z + QExiz  +  c,  O"^  •dsH  =  (A/)^  •  (Savg{z)  ■  Ex(z,  t)  (2.39) 

The  second  curl  equation  of  Maxwell  can  now  be  replaced  by  an  approxiinate  equation 
\Hy{z  -  f ,  0  - H,{z  +  f ,  t)]  ■  A/  « |{  [s^iz)Ex(z,  t)]iArf  ]  +  (AI)^-  o^,(z)  •  E^(z,  t)  (2.40) 
vdiich  can  be  simplified  (because  of  media  stationarity)  to 

Hy{z  -  f ,  0  -  Hy{z  +  f ,  0  «  Za^iz)  ■  A/  •  |{£.(z,  /)}  +  A/  -  a^^(z)  -  Ex{z,  t)  (2.41) 


The  approximate  equation  can  be  re-written  as 

The  above  equation  is  integrated  with  req)ect  to  time  to  get  an  approximate  equation  for 
the  electric  field  at  the  present  time  t 

A  .gmilar  equation  can  be  written  for  the  electric  field  at  the  same  spatial  location 
but  at  an  earlier  time  t-At 

Subtracting  the  two  equations  we  get 

Ex(z,t)^Ex{z,t-At)+  (2.45) 

1 


Al■Za^g(z) 


;l  -  f .  - L  + f  -  si 
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Assuming  that  the  interval  At  is  small  enough  such  that  the  magnetic  fields  can  be 
assumed  approximately  constant  within  At  and  equal  to  tihe  value  at  the  center  of  the  time 
interval  At 


ws  get 


H(z  +  f,x)'^ 

^H(z  +  f,t-f)  for  t-At<x<t 

(2.46) 

H(z-f,x)f 

a//(7-f,/-f)  for  t-At<x<t 

(2.47) 

Exiz,f)!^Exiz,t-Af)+ 

(2.48) 

Ai  r  ET  /_  A/  ^  tj  ,  A/  ^  A^^^  P 

+  ,)J-  I 

The  above  equation  can  be  written  in  a  someA\fiat  differoit  form 

Ex(z,  t)  «  Ex(z,  t  -  At)+ 

At  V  TJ  A/  .  A/X  tr  ^  *  A/xl  ^^aavg(z)r  1 


(2.49) 


The  term  in  the  brackets  is  recognized  as  the  average  value  of  the  electric  field 
within  the  interval  (t-At,  t).  The  above  equation  can  be  written  in  a  more  con:q)act  form 


using  the  identity 


1  _  1  1  _  1  1  _  1  /H£±  =  v  Z  ^ 

e  eoSr  yti^^VSoe.  «  "s. 


(2.50) 


wfiere  the  free  space  velocity  of  propagation  is  denoted  v^,  the  intrinsic  impedance  of  free 
space  is  denoted  (Zo  »  317Q)  and  s,  denotes  relative  permitivity.  Introducing  the  grid 
"propagation"  velocity 

(2-51) 

we  can  write  the  approximate  equation  for  the  electric  field  "updates"  as 

Ex(z,f)«Exiz,t-At)+  (2.52) 


Vq  1 

^grid  Srxzvgiz^ 


Hyiz  -  f ,  t  -  f )  -  Hyiz  +  f ,  t  -  f )  -  A/  •  Qa.g(z)  Ex(z,  x)ck^ 
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or 

E,{z,t)«E^(z,t-M)+  (2.53) 

Vgrid  8  r  mgiz)  ~  ^  ~  T)  ~  +  Y5^“'2')“^^'  <^avgi^Ex,avg(t^)(z,  ?)] 

The  average  of  the  electric  field  between  t-At  and  t,  denoted  can  not  be 

calculated  exactly  because  the  exact  ten^oral  variation  of  the  electric  field  is  not  known  a 
priori.  However,  assuming  a  linear  variation  within  At  interval,  the  average  field  is  the 


average  of  the  values 


at 


(2.54) 


!  E _ ..t)(z)  into  the  update  equation  we  get 


Substituting  the  time-average  E^^^,)( 

Ex(z,  t)  « Ex{z,  t -  Ar)+  (2.55) 

'7  ^0  1  nr  /_  A/  >  AZ-v  IJ  ^  f  Af\  at  (_\  A/) +£j;(z, /)"j 

-  y , y)  -  +  y , «  -  y)  -  A/  •  - ^ - J 

which  gives  the  final  equation  for  the  electric  field  updates  as 

j  1  Vq  Zq  '  A/  •  Ciavg(z) 

(2.56) 

J  1  ^0  ^0  '  '  ^crvg\^) 

2  £  r^ovg 

Vq  1 

^grid  er^vg(z)  f  rr  /_  AZ  A^^  zi  ^  Ar^■l 

— -  y,  <- y) -/?,(.  + y.  <- y)J 

2  '^grid  Sr,avg(z) 

The  equation  simplifies  greatly  for  the  case  of  non-conductive  media  (ct=0) 

EAz,  0  »  £.(z,  /  -  A()  +  Zo  -  f ,  /  -  f)  -  /(,(z  +  f ,  (  -  f )]  (2. 57) 

hr  the  case  of  firee-space  (s  =1),  the  equation  simplifies  further 

t-At)+  Zo  -  f ,  t  -  f )  -  Hy(z  +  f ,  ?  -  f )]  (2. 58) 


Ex(z,t)  «iEx 


19 


These  electric  field  "update"  equations  also  show  that  the  electric  field  E  and  the 
magnetic  field  H  should  be  evaluated  at  points  shiJfted  spatially  by  AI/2,  and  at  instants 
separated  in  time  by  At/2. 

B.  TRANSPARENT  GRID  TERMINATION  IN  1-D 

1.  1-D  Grid 

The  electric  and  magnetic  field  update  equations  have  been  derived  using  local 
coordinate  systems,  with  origins  at  the  centers  of  the  magnetic  and  electric  contours, 
respectively.  These  equations  now  need  to  be  "converted"  to  a  global  coordinate  system, 
that  is  to  the  grid  of  equi-distant  sampling  points  along  the  z-axis  (for  1-D).  We  will 
assume  that  our  domain  ("grid")  is  a  line  segment  of  length  L.  The  fields  are  sanq)led 
using  a  spatial  step  Al^ITN^.  The  electric  and  magnetic  field  sanq)le  locations  are 
"interleaved",  as  shown  below  for  N^=5.  The  first  and  the  last  spatial  sampling  points 
form  the  grid  edges.  The  spatial  edge  samples  can  be  either  electric  field  samples  or 
magnetic  field  samples.  Although  the  selection  of  the  field  for  the  grid  edges  makes  no 
difiFerence  in  principle,  we  will  in  general  use  electric  field  samples  for  the  grid  edges. 

f . O . • . O . • . O . • . O . • . o . • 

A1 

< . > 

Figure  2.5.  1-D  E-field  (black)  and  H-field  (white)  nodes. 
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We  would  evaluate  the  fields  within  the  grid  first,  since  a  known  incident  field 


propagates  from  the  center.  Electric  field  update  equations  can  now  be  "converted"  to 

electric  field  grid  equations  by  replacing  the  variables  z  and  t  with  the  grid  coordinates  of 

the  electric  field  spatial  and  tenq)oral  sairq)]ing  points 

z  ->  M/,  ^  =  0,  \...Nz  and  t  ->•  nAt,  «  =  0, 

The  electric  field  update  equation  for  non-conductive  media  is 

Ex(kAl,  nAt)  «  Ex(kAl,  nAt  -  At)  (2.59) 

The  notation  can  be  sinq)lified  further  by  omitting  the  common  A1  and  At  terms,  and  using 


a  superscript  for  the  index  of  the  ten^oral  san^ling  point 


E",(k')«ETKk)+Z^(j£) 


1 


^grid-f  Sr,avg(A) 


H7\k-\)-H7\k+\) 


(2.60) 


The  electric  field  grid  equation  for  conductive  media,  using  the  same  notation  as  above,  is 


E”xik) 


1  \(  Vo  ^0  '  Al  •  Oavg(]^ 


'rfivg 


{k) 


1  I  if  vq  '\Zo-Al-o„,g(k) 
ly^grid)  Zr^gik) 


ET\kyr 


(2.61) 


'ygridJ  Zr^gjk) 


1  if  Vq  ^  Zq  •  A/  •  Cavg(k) 
^'^'yyv^ridJ  p_ 


ff;"(k-^)-ffp(k+h 


2  \VgridJ  Zr,avgiE) 

Finally,  the  free-space  electric  field  grid  equation  is 


£;(i)»£r‘W+Zo(^) 


(2.62) 


Similarly,  magnetic  field  update  equations  can  be  converted  to  magnetic  field  grid 
equations  by  replacing  the  variables  z  and  t  with  the  grid  coordinates  of  the  magnetic  field 
spatial  and  ten^oral  sampling  points 
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z  (^+ |■)A/,  A:  =  0,  \...Nz  -  1  and  t^(n->r |)A?,  «  =  0,  L-.A'/ -  1 
The  magnetic  field  grid  equation,  using  the  same  notation  as  for  the  electric  field  grid 

equation,  is 

Hf\k  +  i)  »  Hf\k + i)  -  Ko  fe) - TT-d-Eifk  + 1)  -  E"m  (2.63) 

or,  for  non-magnetic  media 

Hf\k+\)«Hf'(k+\)-Y^{^[E-.(k+l)-E"m  (2.64) 

The  electric  and  magnetic  field  grid  equations  have  the  following  general  form 

Er  =  CsiEf  +  CEiVHf  (2.65) 

=  CmHf  +  CHT^Ef  (2.66) 

where  C's  are  constants  (real  numbers)  that  depend  on  the  media  properties  and  the 

velocity  ratio  v^/Vg^^ ,  and  "del"  operator  represents  the  spatial  derivative  (gradient).  This 

general  form  of  the  grid  equations  can  be  interpreted  as  follows;  "the  new  value  of  EJH 

field  at  a  grid  node  is  equal  to  the  weighted  sum  of  the  old  value  of  the  E/H  field  at  the 

same  node  and  the  spatial  variation  of  the  old  H/E  field  between  the  two 

nearest-neighbor  nodes. " 

2.  Grid  Termination 

The  grid  equations  derived  above  are  valid  for  all  the  nodes  except  for  the  nodes 
on  the  grid  edges.  The  reason  is  that  the  grid  edge  nodes  have  only  one  neighbor  node 
instead  of  two  like  any  node  that  is  not  on  the  grid  edge.  The  edge  nodes  with  a  single 
"neighbor"  thus  need  to  have  equations  different  than  the  equations  we  have  derived  for 
the  "non-edge"  grid  nodes  with  two  "neighbors".  TGT  requiremait  is  conceptually  very 
simple  and  straightforward  to  m5)lement  in  1-D.  To  that  effect  we  will  use  the  concept  of 
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a  multiport  (a  two-port  in  1-D).  It  does  not  matter,  in  principle,  wliether  the  edge  node  is 
an  E  field  or  an  H  field  node.  We  will  assume  that  the  edge  node  for  the  two-port  model 
is  an  E  field  node.  The  input  port  of  the  two-port  will  be  the  nearest  node  "of  the  same 
kind".  Since  we  have  assumed  an  E  field  node  for  the  output  port,  the  input  port  will  be 
the  nearest  E  field  node  inside  the  grid.  We  could  have  also  selected  an  E  field  node  as  the 
ou^ut  port  and  the  nearest  H-field  node  as  the  input  port.  However,  the  selection  of  the 
same  kind  of  node  for  both  input  and  output  ports  has  the  advantage  that  the  TGT  resuhs 
obtained  in  this  manner  can  be  also  used  to  solve  the  wave  equation  (a  second  order 
partial  differential  equation)  that  has  only  one  field  as  the  variable  and  the  grid  with  only 
one  kind  of  nodes.  The  figure  below  shows  the  two  ports  for  the  two  (1-D)  grid  edge 
nodes. 


out 


h(t) 


in 


-o 


o 


m 


o 


h(t) 


out 


grid  edge 


grid  edge 


Figure  2.6.  Modeling  1-D  Grid  Terminations  as  Two-ports. 

The  fields  at  the  output  ports  E(0,t)  and  E(L,t)  can  be  expressed  as  convolutions 
of  the  fields  at  the  input  ports  E(Al,t)  and  E(L-Al,t)  and  the  two  port  inqiulse  response 
h(t).  (The  impulse  responses  for  the  two-ports  on  the  left  and  on  the  right  edges  are 
identical,  by  symmetry). 

£,(0, 0  =  i)  *  hit)  =  fo  £,(A/,  x)hit  -  x)di  (2.67) 

E,(L,  t)  =  E^iL  -  A/,  0  *  hit)  =  j;  E,iL  -  A/,  x)hit  -  x)dx  (2.68) 
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The  discretized  forms  of  these  equations,  involving  saroples  of  the  jfields  at  t=nAt,  are  the 
discretized  convolutions 

=  (2.69) 

E"(N,)  =  25  EiiN,  -  1)  •  h"-P  (2.70) 

The  convolution  equations  ejq)ress  the  fields  at  the  edges  as  weighted  sums  of  the 

time  histories  of  the  fields  "just  inside"  the  grid.  In  that  respect  they  are  equations  of  the 

same  type  as  the  equations  for  the  non-edge  grid  nodes,  except  that  they  involve,  in 

general,  summations  with  more  terms.  However,  the  inqrulse  response  h(t),  as  will  be 

shown,  is  a  rapidly  converging  (to  zero)  fimction  which  reduces  the  number  of  relevant 

terms  in  the  convolution  sum.  The  impulse  response  h(t)  (actually  its  "san:5)led"  form  h“) 

needs  to  be  determined  only  once,  for  a  selected  grid  velocity  Vgrij=Al/At.  The  issue 

remains  how  to  determine  the  hiq)ulse  response?  The  discretized  iaq)ulse  response  h“  will 

be  determined  using  the  discrete  equivalent  of  the  Dirac  deha  fimction  which  we  will 

denote  as  5“  Since  there  are  two  grid  edges  and  two  input  ports  (Fig.  2.6)  one  input  port 

will  be  set  to  zero  and  5”  will  be  applied  to  the  other  iD|)ut  port.  Since,  in  1-D,  there  is 

only  one  h”  to  determine  we  will  set  E,j(Al,t)  to  zero  and  apply  the  Dirac  deha  fimction  as 

E,(L-Al,t) 

E"(N,  -  1)  =  5"  (2.71) 

i;;(l)  =  0  (2.72) 

This  is  depicted  in  Figure  2.7.  Note  that  the  mpulse  response  h°  is  observed  at  z  = 

L  and  that  the  grid  extends,  theoretically,  to  infinity  past  the  observation  point  z  =  L. 

Since  5“  =  0  for  n  >  0,  the  grid  to  the  left  of  the  (L  -  Al)  node  is  effectively  isolated  form 
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the  observation  point  at  L.  This  means  that,  to  find  h",  we  may  just  consider  the  grid  to 
the  right  of  the  L-Al  node,  as  shown  below. 

^Al,t)=0  E(L.Al,t)=5(t) 

^  h(iO 

# . O . • . O . • . O . # . O . • . O . • . o . t . O . • . O’  • 

A  A 

gridec^  gridecfee 

E^I^Al,t)=5(t) 

•  © • . ©--•-© . - . • n . • 

A  w  w 

gride^ 

Figure  2.7.  Determining  the  Irq)ulse  Response. 

The  impulse  response  needs  to  be  obtained  as  if  the  grid  were  not  terminated  at  all 
to  the  right  of  the  observation  point  at  z  =  L.  If  the  grid  were  terminated,  the 
"reflectionless"  termination  condition  would  be  needed  and  this  is  exactly  vdiat  we  do  not 
have  and  are  trying  to  find.  A  grid  extending  to  infinity  does  not  require  a  termination 
and  thus  "avoids"  the  termination  problem  Although  it  is  not  pos^le,  in  practice,  to 
extend  the  grid  to  infinity,  the  grid  can  be  made  large  enou^  such  that  any  "reflections" 
off  the  new  grid  edge  would  have  arrived  after  the  impulse  rehouse  has  "converged"  to  a 
small  selected  value  that  we  consider  as  "zero".  The  time  it  takes  the  iirpulse  to 
"propagate"  to  this  new  grid  edge  and  back  to  the  observation  point  is 

Td=^  =  ^^^^  =  2NDAt  (2.73) 

At 
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where  D  is  the  distance  between  the  observation  point  and  the  new  grid's  edge  and  Nj,  is 
the  number  of  nodes  between  the  observation  point  and  the  new  grid's  edge.  The  duration 
of  the  iiqpulse  response  h“  that  will  not  be  "corrupted"  by  the  reflections  will  be  2Np. 
Finally,  it  is  assumed  that  the  space  outside  the  grid  is  free-space,  and  the  grid  equations 
for  the  electric  and  magnetic  fields  will  thus  be  the  free-space  equations 

£;(i) =£:;-'(*) +Zo(^gj)  (2-74) 

<275) 

The  order  of  the  terms  in  the  brackets  for  the  magnetic  field  equation  has  been 

reversed,  such  that  the  electric  and  magnetic  field  grid  equations  will  have  the  same  sign 

(+)  in  front  of  the  brackets.  The  above  equations  can  be  also  written  in  a  different  form, 

to  reduce  the  number  of  multipUcation's  that  would  need  to  be  done  at  each  time  step. 

Multiplying  the  magnetic  field  equation  by  Z^,  we  get 

ZoHf\k+\)  =  ZoH7^\k+\)  +  (^)  1)]  (2.76) 

An  auxiliary  variable  hy  may  be  introduced 

hy  =  ZoHy  (2.77) 

and  the  electric  and  magnetic  field  grid  equations  can  be  written  using  h^,  instead  of  Hy 

E'i(k)=ETKk)  +  {j^  hf\k-\)-hfHk+^)  (2.78) 

hf'  (4  +  0  =  h7\k + i)  +  (^)  [£;(i)  -  £«*  +  01  (2.79) 

The  equations  for  the  electric  field  (E^  )  and  for  the  magnetic  field  multiplied  by 
(ZjHy)  have  identical  forms,  involving  only  one  parameter,  the  velocity  ratio  The 

boundary  inopulse  response  h“,  to  be  obtained  using  these  equations,  will  thus  be  valid  only 
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for  a  selected  velocity  ratio  or,  equivalently  for  the  selected  grid  velocity  Vgrij=Al/At  The 
boundary  in^ulse  response  h”,  obtained  using  the  equations  for  V(/Vjj,jj=l,  is  shown  below. 
Note  that  the  plotting  program  interpolates  linearly  the  values  between  the  san:q)ling  points 
nAt. 


Time  Step 

Figure  2.8.  Boundary  Inq)ulse  Response. 

The  boundary  iiiq)ulse  response  for  1-D  is  very  sinq)le 

h"  =  51  (2.80) 

The  in:q)ulse  response  is  non-zero  only  for  n  =  1  where  it  has  the  value  of  1,  whidi  is  a 

delta  inq)ulse  delayed  by  At  and  may  also  be  written  as 

/2(0  =  5(/-A0  (2.81) 

We  can  now  write  the  grid  equations  for  1-D  grid  edges  as  follows 

^?(o)=i?:s^?(i)-A”-^=£:r'(i)-/2'  =i:r'(i)  (2-82) 

E%N,)  =  EI(N,  -  1)  •  =  ETKN,  -  1)  •  hi  =  ETKN,  - 1)  (2.83) 
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The  above  equations  state  that  the  fields  at  the  grid  edges  are  updated  by  siroply 
taking  the  previous  values  of  their  nearest  neighbors  inside  the  grid.  Although  the 
boundary  iropulse  responses  are  not  so  sinaple  for  2-D  and  3-D,  they  can  be  determined 
using  essentially  the  same  procedure  as  shown  for  1-D. 

C.  1-D  TGT  RESULTS 

We  apply  the  geometrical  model  shown  in  Figure  2.6.  The  nodes  on  the  edges  are 
E  field  nodes.  The  source  is  the  electric  field  at  the  center  node  and  the  source  waveform 
is  a  unit  an^litude  delta  pulse .  This  represents,  in  1-D,  a  uniform  plane  wave  propagating 
fi’om  the  center  into  +z  and  -z  directions.  By  applying  the  "standard"  FD-TD  equation  for 
the  nodes  inside  the  grid  and  the  TGT  equations  for  the  grid  edges  we  obtain  the  power 
within  the  grid  as  a  fimction  of  time  as  shown  below.  It  is  clear  that  the  power  within  the 
grid  is  constant  until  the  wave  has  "left"  the  grid  vshen  it  fells  to  zero  abruptly.  The  1-D 
grid  with  TGT  termination  thus  behaves  like  an  "ideal”  grid. 
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Figure  2.9.  Residual  Power  for  TGT. 

We  will  next  con^are  this  result  with  that  for  an  "infinite"  grid.  Note  that  the 
power  is  calculated  within  the  same  grid  region  as  in  the  previous  Jfigure. 


Figure  2.10.  Residual  Power  for  an  "Infinite"  Grid. 
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in.  ANALYSIS  OF  TGT  FOR  2-D  FDTD 


Figure  3.1.  TE^  Field. 

The  direction  of  propagation  is  indicated  by  the  propagation  velocity  vector  v  . 
The  unit  vector  in  the  direction  of  propagation  and  the  electric  and  magnetic  field  vectors 

form  a  triplet  of  mutually  orthogonal  vectors.  The  "dual"  TM^  field  has  the  conq)onents 

H^H^andE, 

^ m(x,y, t)  =  Hx(y,y, t)x  -\-Hy{x,y, t)y  (3.4) 

E  TM(x,y,  t)  =  Ez(x,y,  i)z  (3.5) 

as  shown  below. 

H 

Hx . - . - . -••>  X 

Figure  3.2.  TM,  Field. 
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Our  objective  is  to  determine  the  discretized  forms  of  Maxwell's  curl  equations  in 
the  integral  form  for  TE^  fields 

^Ex(x,y, iyx  +Ey{x,y, t) •  dh  =  V-{^,y)Hz{x,y, t)~z  •  dsE  ■  (3.6) 

\Hz(x,y,  tyz~^»dlH  =  s(^,  t)x  +Ey{x,y,  t)y  ]  •  ^3^s//|+  (3.7) 

0^  +Ey(x,y,  t)y  ]  •  dsn 

and  for  TM^  fields 

\^Ez{x,y, t)z  ]•£//£  =  -^1  \^(y,y\Hx(x,y,  f)x  +Hy{x,y, t)y  ]  •  •  (3.8) 

[//;,(x,3;, 0^ +Hyix,y,  t)y  ]  •  =  (3.9) 

^|ILh  '\*^n-  +Jf^^  o(x,y'^Ez(x,y, tfz  '^•dsn 

Note  that  the  the  second  curl  equations  do  not  have  the  source  current  terms,  since 
it  has  been  assumed  that  there  were  no  somce  currents  in  the  domain  of  interest.  The 
contours  of  integration  for  the  electric  and  the  magnetic  field  circulations  are,  in  general, 
different  and  are  thus  labeled  Cg  for  an  electric  field  contour  and  for  a  magnetic  field 
contour.  Similarly,  the  surfeces  associated  with  the  contours  are  labeled  Sg  for  the 
magnetic  flux  and  Sjj  for  the  electric  flux.  The  electric  and  magnetic  field  contours  for  a 
TE^  field  are  shown  below. 
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Ex(x,jMly/2,t) 

, . — ;^J.' . 

Ey|;x-dx/2,y,0  ^  ^Ey(x+dx/2,y,t) 

. 

. . L, . .  . 

Ex(x,y“d)^^) 


Figure  3.3.  Contours  for  TE^  Electric  and  Magnetic  Field  Circulations  and  Fluxes. 

The  electric  field  contours  C^  and  the  magnetic  field  contours  C^  are  square  (A1  by 
Al)  contours,  but  in  orthogonal  planes.  The  Cg  contours  are  in  the  xy-plane,  while  the  Cjj 
contours  are  in  the  yz-plane.  The  contours  can  be  compared  to  links  of  a  (2-D)  chain 
fence,  evoking  the  idea  of  electric  and  magnetic  field  linkage  in  two  orthogonal  directions. 
The  "dual"  Cg  and  Cjj  contours  for  a  TM^  field  are  diown  below 


Hx(x,y-dy/2,t) 


Figure  3.4.  Contours  for  TM^  Electric  and  Magnetic  Field  Circulations  and  Fluxes. 


The  electric  and  magnetic  contours  and  their  associated  surfaces  will  be  used  to  discretize 
Maxwell's  curl  equations. 


2.  Derivation  of  TE^  Magnetic  Field  Update  Equation 

The  discrete  equivalent  of  the  TE^  electric  field  circulation  will  be  determined  next. 
A  contour  Cg  is  shown  below 


Ex(x,y+dy/2,t) 


Ey(x-dx/2,y,t) 


1  ^  I  Ey(x+dx/2,y,t) 

I  Hz(x,y,t)  I 


Ex(x,y-dy/2,t) 


F^ure  3.5.  A  Contour  for  TE^  Electric  Field  Circulation  and  Magnetic  Flux  Calculations. 

The  center  of  the  surface  is  assumed  to  have  the  coordinates  (x,  y).  We  will 
assume  a  muform  grid  with  Ax=Ay=Al.  A  local  coordinate  system  v|/)  will  be 
estabhshed,  with  the  origin  at  the  center  of  the  contour  such  that  any  point  (x'jy*)  within  or 

on  the  contour  Cg  can  be  specified  by  its  local  coordinates  ^  and  v|; 
x'=x  +  ^  y=y  +  v}/  where  and 

The  local  coordinates  will  be  used  in  evaluation  of  the  line  and  surface  integrals 
that  constitute  the  integral  forms  of  Maxwell's  equations.  The  circulation  of  the  electric 
field  around  the  electric  field  contour  C^  ,  assuming  counter-clockwise  reference  direction 
such  that  the  normal  to  the  surface  is  in  the  direction  of  the  magnetic  field  H(x,y,tX 
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can  be  evaluated  approximately,  assuming  that  the  electric  fields  are  constant  over  the 
length  A1 ,  along  the  edges  of  the  (A1  by  Al)  square  contour 

t)x  +Ey{x',y', t)y  ]•£//£«  (3. 10) 

[E,(x,y  -  f ,  0  +Ey{x  +  f  t)  -  E,ix,y + f ,  0  -  Ey{x  -  f,y,  0]a/ 

The  magnetic  flux  through  the  surface  can  be  calculated  using  the  local  coordinates 

(3.11) 

+ %.y + V,  O’?  • 

There  are  infinitely  many  ways  to  model  the  variation  of  the  magnetic  flux  density 

with  the  local  coordinates  ^  and  \|/  within  the  contour.  We  first  need  to  postulate  a 

certain  variation  of  the  magnetic  field  with  ^  and  v|/  such  that  the  integral  can  be  evaluated 

over  Sg.  The  sinqjlest  model  assumes  that  the  contour  width  Al  is  small  enough  such  that 

the  magnetic  field  H(x+^,y+\|/,t)  within  the  contour  may  be  assumed  constant  and  equal  to 

the  value  at  the  contour's  center  H(x,y,t). 

Hz(x  +  ^,y  +  \V,t)>=iH^(x,y,t)  (3.12) 

This  is  equivalent  to  using  a  piece-wise  constant  or  2-D  "pulse"  expansion 

approximation  of  the  actual  magnetic  field  variation  with  the  z-coordinate.  The  above 

assumption  allows  that  magnetic  field  ,  although  constant  within  a  contour,  can  change 

firom  one  contour  to  an  other.  This  yields  an  approximate  e?qpression  for  the  for  the 

magnetic  flux 

IL  yWzix' y ,  iyz  .  «  Hz(x,y,  t)  •  0  0  p(x + l,y  +  (3.13) 

Se  2  2 

The  integral  of  the  permeability  can  be  re-written  in  the  following  manner 
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0  0  0  0  [x(x + ^,y  +  \\f)d^d\\f  (3. 14) 

2  2  1_  (A/)  2  2 

The  term  in  the  brackets  is  recognized  as  the  average  permeability  within  the  contour 

Ce 

=  (3.15) 

(A/)  2  2 

The  approximate  expression  for  the  magnetic  flux  through  can  now  be  written  using  the 

average  permeability  as 

J0  ii(x'y)Hz(x'y,  ifz  •dsE~  {Mf  \i.avg{xyH^(x,y,  i)  (3. 16) 

Again,  this  approximate  e?q)ression  resulted  fi'om  the  piece-wise  constant 

approximation  of  the  magnetic  field  with  respect  to  x  and  y-coordinates.  The  main 
advantage  of  the  piece-wise  constant  expansion  employed  above  is  its  simphcity.  Better 
accuracy  can  be  achieved  by  using  more  involved  models  for  the  field  variation  with  x  and 
y,  but  at  the  expense  of  increasing  the  con^utational  time.  The  first  ciul  equation  of 

Maxwell  can  now  be  rqrlaced  by  an  approximate  equation 

[Ex{x,y-^,t)+Eyix-k-^,y,t)-E:,(x,y  +  f,t)-Ey{x-^,y,t)\-M «  (3.17) 

which  can  be  sinq>lified  as  below  because  of  media  stationarity 

Ex(x,y-^,t)+Ey{x  +  ^,y,t)-Ex{x,y  +  ^,t)-Ey{x-^,yJ)«  (3.18) 

-jf{Hz(x,y,  0}  •  A/  -  fiavg(x,y)} 

This  equation  can  also  be  re-written  as 

|{^.(x,:^,0}  (3-19) 

A/-^^(xy)[^''^^’-^  -  f,  t)  +Ey{x + f  ,y,  /)  -  Ex{x,y + f ,  t)  -  Ey(x  -  ^,y,  r)] 
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The  time;  derivative  operator  in  the  above  equation  is  typically  "replaced"  by  the 
finite  difference  approximation.  Ho\vever,  just  like  in  1-D,  we  present  an  alternate 
approach,  such  that  the  approxunation  of  the  field  time  variation  is  shown  to  be  analogous 
to  the  approximations  already  introduced  for  the  field  spatial  variation.  The  above 
equation  is  integrated  with  respect  to  time  to  get  an  approximate  equation  for  the 


magnetic  field  at  the  present  time  t 


(3.20) 


[ j'  E,(x,y  -f.x)dx-  J'  E,(x,y + f .  t)A  +|'  £,(x  +  f,y,  iyk  -  J'  £,(x  -  f t)*] 

A  similar  integral  equation  can  be  written  fijr  the  magnetic  field  at  the  same  spatial 


location  but  at  an  earlier  time  t-At 


(3.21) 


Subtracting  the  two  equations  we  get 


H,(x,y,  t)  «  H,(x,y,  t- At)  - 


(3.22) 


[t^E,<x,y-^.zyh-l_^E,{x,y+^.z)dxyl_^E,fx+i,y.zyk-^^^E,(x-!^,y.z)ch] 

The  integrals  on  the  right-hand  side  can  not  be  evaluated  exactly,  because  the 
exact  temporal  variation  of  the  electric  fields  within  the  At  interval  prior  to  t  is  generally 
not  known.  However,  the  integrals  can  be  evaluated  approximate^  by  assuming  a  certain 


variation  of  the  electric  field  with  the  temporal  variable  t.  The  amplest  approach. 


consistent  with  the  assumptions  made  for  the  field  spatial  variation,  would  be  to  assume 
that  the  interval  At  is  small  enough  such  that  the  electric  field  can  be  assumed 
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approximately  constant  within  At  and  equal  to  the  value  at  the  center  of  the  time  interval 
(t-At,  t) 


-  f .  i)  “  E,(x,y  -  f ,  /  -  f ) 

for 

t-At<x<t 

(3.23) 

E.(x,y+f,x)«E.(x,y  +  f,t-f) 

for 

t-At<x<t 

(3.24) 

1 

1 

1 

for 

t-At<x<t 

(3.25) 

E,(x+f,y,t)«Ey(x+f,y,t-f) 

for 

t-At<x<t 

(3.26) 

The  above  represents  a  piece-wise  constant  approximation  of  the  electric  field  variation 
with  respect  to  the  ten:q)oral  variable  t.  The  approximate  expresaon  for  the  magnetic  field 


at  location  (x,y)  and  at  the  time  t  now  becomes 

(3.27) 

The  above  equation  can  be  written  in  a  somewhat  different  form,  using  the  identity 

2. _ 1 _ 1 _ _ 1 _ L  —  1 _ 1__L  —  Vn-i _ =  Vn  Tn—  (3  28t 

P-r  JJx^  °Zo  1^-  °  °  1^- 


vsdiere  the  free  space  velocity  of  propagation  is  denoted  Vj,  the  intrinac  mq>edance  of  free 
space  is  denoted  (Zo  «  377fl),  and  p,  denotes  relative  permittivity.  Introducing  the 
grid  "propagation"  velocity 


_  A/ 


(3.29) 


we  can  write  the  approximate  "update"  equation  for  the  magnetic  field  as 

i)  ^ H,(x.y, t-At)-  « 


(3.30) 


[E:,(x,y-f,t-j)-  E:,{x,y  +  f ,  /  -  f )  +  £^(x  +  f  ,y,  /  -  f )  -  Ey{x  -  f,y,  t  -  f  ] 


The  equation  suiq)lifies  for  the  case  of  non-magnetic  media  (Pr“0 


^fz(x,y,  t)  »  Hz(x,y,t-At)  -  x 


(3.31) 


[&(x,:v  -  f ,  ( -  f )  -  E.(x.y + f  ■ '  -  f )  + f  .T-. '  -  f  > -Ey(x  -  f  ,3>. '  -  f)] 
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The  last  two  equations  show  that  the  magnetic  field  H  and  the  electric  field  E  are 
evaluated  at  points  shifted  spatially  by  Al/2,  and  at  instants  separated  in  time  by  At/2,  just 
like  in  1-D.  The  relationship  between  spatial  and  tenq)oral  "sanq)les"  of  the  electric  and 
the  magnetic  fields  is  thus  the  same  and  the  samples  are  shifted  with  respect  to  each  other 
by  one-half  of  the  sanq)ling  interval.  The  approximate  equations  for  the  electric  field 
"updates"  can  be  derived  using  the  same  procedure  shown  for  the  magnetic  field  update 
equations. 


3.  Derivation  of  TE^  Electric  Field  Update  Equations 

We  will  start  with  the  Maxwell's  curl  equation  for  the  circulation  of  the  magnetic 


field 


Hz{x,y, t)~z  •  dl a  =  t)x  +Ey(x,y, /) 7  ]  •  ^sh 

(s(x,y)  Ez(x,y, t)x  +Ey(x,y, t)y  ] • 


+  (3.32) 


The  discrete  equivalent  of  the  magnetic  field  circulation  needs  to  be  determined 
for  two  sanq)le  contours  a  contour  in  a  plane  parallel  to  the  yz-coordinate  plane,  and  a 
contour  in  a  plane  parallel  to  xz-coordinate  plane. 


Hz(x-A  l/2,y,t) 


.>x 


Hz(x  +  Al/2,y,t) 


Figure  3.6.  A  Contour  for  Updating  E^. 
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A  local  coordinate  system  (^,  Q  will  be  used,  with  the  origin  at  the  center  of  the 
contour.  Any  point  (x',z')  within  or  on  the  contour  Cjj  can  be  specified  by  its  local 


coordinates 


x^=^x  +  %  z'  =  z  +  l^  where  -  —  <^<4-;^  and 


2 


The  circulation  of  the  magnetic  field  around  ,  assuming  counter-clockwise  reference 

direction  such  that  the  normal  to  the  surface  Sjj  is  in  the  direction  of  the  electric  field 
Ey{x,y,  t),  is  therefore  given  exactly  by  the  following  simple  expression  (because  there  is 
no  magnetic  field  variation  with  the  z-coordinate) 

I  H,{x',y',  fit  =  0  H,{x -  f ,y,  fft  •  -  0  +  ^,y,  fit  •  (3.33) 


f  H^(x',y',  t)t*dlH  =  [h,{x  -  f,y,  f)  -  H,(x  +  f,y,  t)]  ■  A/  (3.34) 

The  electric  flux  through  the  surfece  can  only  be  evaluated  approximately,  since 
the  exact  way  that  the  electric  flux  denrity  Dy(x,y,  t)  =  zix,y')Ey{x,y,  t)  varies  with  x  and  y 
is  not  known  a  priori.  First  we  need  to  evaluate  the  flux  integrals 


0^  s(x',y,t)Ey(x',y',  fiy  *dsH  = 


(3.35) 


0^  o(x',y')Ey(x',y',  tiy  •dsH  = 

0  0  cy(x+4,y +  vi;)£;c(x+^,3^  +  v|/,  fiy  •  'yd^^dC, 


3.36) 


Since  the  integrands  are  not  fimctions  of  the  local  z-coordinate  C  3nd  the  local 


coordinate  h/=0,  the  surfece  integrals  are  sinq)lified  to  line  integrals 

0^  8(x + ^,y)Ey(x  +  ^,y,  fiy  •dsH==M  0  s(x + %,y)Ey{x  +  l,y,  t)^  (3.37) 
ff  (j(x + 'i,y)Ey{x  +  ffy  •dsM  =  M  0  cy(x + l,y)Ey(x + ^,y,  t)d^  (3.38) 

JJi>H  •'-j 
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Next,  we  assume  that  the  contour  side  A1  is  small  enough  such  that  the  electric 
field  E(x+^,y,t)  within  the  contour  naay  be  assumed  constant  and  equal  to  the  value  at  the 
contour's  center  E(x,y,t).  This  assumption  yields  approximate  expressions  for  the  flux 


integrals 

ff ^  s(x + ^,y)Ey(x  +  ^,y,  t)'y  Ey(x,y,  t)  ■  H  (3-39) 

Sh  2 

(3{x  +  l,y)Ey{x  +  '^„y,  t)y  •dsH»  M  Ey(x,yJ)  -0  Q(x  +  '^,y)d^  (3.40) 

S//  2 

or,  using  the  permittivity  and  conductivity  averaged  in  the  x-direction 

ff  s(x  +  ^,>’)£:/x+iy,07*3/,«(A0"  -Ey{x,y,t)-z^g,(x,y)  (3.41) 


IT  q(x  +  i^^y)Ey(x  +  ^,3^,  i)y  •  dsn  «  (A^  •Eyix^y^  t)  •  cjav^xO^yy)  (3.42) 

JJS/i 

The  second  curl  equation  of  Maxwell  can  now  he  replaced  by  an  approximate  equation 

[H,ix  -  f  ,y,  0  -  H,(x  +  f ,  j,  O]  -  A/  «  (3.43) 

^  { [Za.gx(x,y)Eyix,y,  }  +  {Alf  ■  cJ«^(x,y)  •  Ey{x,y,  t) 

vviiich  can  be  simplified  to 

H,(x-f,y,  t)-H,(x+f,y,t)  «  (3.44) 

£aygjr(^5y)  *  A/ •  -^^Ey(x,y,t)^  +  A/ •  <y<nyit(x,_y)  •  Ey(x,y,f) 


The  approximate  equation  can  now  be  re-written  as 

;>,()>»  (3-45) 

—L—[H,(x  -  f.y,  0  -H.(xy  f,y.  r)]  -  •  E,(x,y.  i) 

Al  ■  Zavgx(x,y)'-  ^  -* 

The  above  equation  is  integrated  with  reject  to  time  to  get  an  approximate  equation  for 


■Ey(x,y,t) 


the  electric  field  at  the  present  time  t 
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A  gitnilar  equation  can  be  written  for  the  electric  field  at  the  same  q)atial  location  but  at  an 
earlier  time  t-At 

Ey(x,y,t-At)  «  A/ •  f  f (3-47) 

^a.g.(x,y)  Jo 


Subtracting  the  two  equations  we  get 

^y(x,y,  0  -  ^y(x,yy  f  -  At)  » 


(3.48) 


which  can  also  be  written  as 


Ey(x,yA 


(3.49) 


Ey(x,y,  t-Ai)  +  i^[  H,(x  -  f  ,7,  t)  A  -  H.(x + f,y,  x)di] 

_gavg.(x,j|p  Ey(x,y,T)ck 

Savg.(x,y)  3>^  ^ 

If  the  interval  At  is  small,  the  magnetic  fields  can  be  assumed  approximately  constant 


within  At  and  equal  to  the  value  at  the  centw  of  the  time  interval  At 


H,(x  +  f,y,x)^H,(x  +  f,y,t-f)  for  t-At<x<t 
H,(x-f,y,T)»H^x-^,y,l-f)  for  t-At<zil 


(3.50) 

(3.51) 


Using  the  above  approximation  we  get 
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(3.52) 


Ey(x,y,t-M)  + 


Ey(x,y,t)« 

{ft(x  -  f,y,  + 


A/  •  Savgx(x,y) 

(x,y) 


^avgxiX^y) 

The  above  equation  can  be  written  in  a  somewhat  dififerent  form 

Ey{x,y,t)K 


(3.53) 


Eyix,y,t-M)  + 


At 


Al  •  Savgx(x,y) 


[h,(x  -  f,y,  ( -  f )  -  ft  (X  +  f  -  f )] 


The  term  in  the  brackets  is  recognized  as  the  average  value  of  the  electric  field  coirq)onent 

Ey  at  (x,y)  within  the  interval  (t-At,  t).  The  above  equation  can  be  written  in  a  more 

compact  form  using  the  identity 
1  1  1 


1 


1 


1 


_ ^ _ ^ _  H£i-  =  v  Z 

e  eoer  V  so  s.  o  <>8. 


(3.54) 


vriiere  the  firee  space  velocity  of  propagation  is  denoted  v^,  the  intrinsic  in:q)edance  of  firee 
space  is  denoted  Z,,  (Zo  «  377Q)  and  s,  denotes  relative  permittivity.  Introducing  the  grid 
"propagation"  velocity 

=  (3.55) 

we  can  write  the  approximate  electric  field  "update"  equation  as 

Ey(x,y,t)!!i  (3.56) 

Ey(x,y,  t-At)+ Zo  Sr.m>?!(x,y)  ~  ^ 

The  average  of  the  electric  field  between  t-At  and  t,  denoted  E,^^,)  can  not  be 


calculated  exactly  because  the  exact  temporal  variation  of  the  electric  field  within  the  At 
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interval  is  not  known  a  priori.  However,  assuming  a  linear  variation  within  At  interval,  the 


average  field  is  the  average  of  the  values  at  t-At  and  t 


.£>,avg(A/')0''5.yj  0  ' 


Ey(x,y,  t-  At)  +Eyix,y,  t) 


(3.57) 


Substituting  the  time-average  ^^,)(x,y,t)  into  the  update  equation  we  get 


Ey{x,y,t)^  (3.58) 

Ey(x,y,  t  -  At)  +  Zo  “  f ^  -  f)  “  +  f ^  “  f )]“ 

Vo  1  _  .  .  Eyix,y,t-M)+Ey{x,y,t) 

which  gives  the  final  equation  for  the  electric  field  updates  as 

J  1  Vq  Eq  -  a/  •  crai^(x,y) 

E^x,y,t,’‘  y  Ey(x,y, I- At}+  (3.59) 

I  _j_  1  "^0  -^0  •  '  ^avgxix^y) 

I'^grid  S;.,<^X(X,>') 

Z  ^0  1 

^grid  Sr^gx(x,y') 

r  TT  Id  f  6t\  fJ  (y-  \  ^  f 

- — — - — 7— "tL^zC^ -  t-  --nzCx  +  Y,y,  t -  y) J 

^  ^  1  Vo  Zo-M-aa^gx(x,y)'- 

2  ^r,avgx(x,y') 

The  equation  sinqrlifies  greatly  for  the  case  of  non-conductive  media  (o=0) 

Ey(x,y,t)x  (3.60) 

Ey(x,y,  t  -  At)  +  Zo  g^^](x,y)  ^  -  f)  “  + f t  -  y)] 

In  the  case  of  fiee-space  (e  =1),  the  equation  sinqrlifies  fiirther 

E/x,y,  t)  »  E,(x,y,  t-At)yZo^[H,ix-  f,y,  t-f)-H,(x+f,y,t-f)]  (3.61) 
The  above  equations  are  the  update  equations  for  the  y-directed  coirqronent  of  a  TE^ 
electric  field.  Similar  equations  will  be  formulated  next  for  the  x-directed  electric  field 


(3.60) 


conq)onent. 
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A 


H2i(x,y+dl/2,t) 


|Ex(x,y,t) 


Hz(x,y-dl/2,t) 


Figure  3.7.  A  Contour  for  Updating  E^. 

A  local  coordinate  system  will  be  used,  with  the  origin  at  the  center  of  the 
contour.  Any  point  within  or  on  the  contour  Cjj  can  be  specified  by  its  local 

coordinates 

y'=y  +  \|/  z'=z  +  ^  where  and  -^<^<-1-^ 

Then,  by  applying  the  same  procedure  drived  for  Ey  field  component,  we  get  the  final 
equation  for  E,,  electric  field  corrponent  as  shown  below 


j  1  Vq  Zq  •  hi  •  O crvgyix,y) 

r.  ^  .N..  2'^  grid  Er/»vgy(x,y)  r  (-,-  x-  f  A  A  I 

Ex(x,y,  t)  « - ^ — 77 - 7 — :r^xKX,y,t-^t)+ 

,  j  vp  Zq  •  A/  •  Oavgy(y,yy 

■‘‘2^  ^rfivgy(x,y) 

Vo  1 

- +  f ,  ( -  f )  -  H.(x,y  -  f ,  <  -  f )] 

1  Vo  Zq  •  A/  •  Oavgy{X,y) 

2  ^grid  ^r,avgy(x,y') 

The  equation  sinqilifies  greatly  for  the  case  of  non-conductive  media  (cr=0) 


(3.62) 
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(3.63) 


E:,(x,y,t)^ 

E,(x,y,  t-M)+  +  f ,  ?  -  f )  -  H,{x,y  -f,t-  f )] 

In  the  case  of  free-space  (8=1),  the  equation  siiiq)Mes  further 

E,(x,y,  t)  «  E,(x,y,  t  -  At)  +Zo:^[H,(x,y  +  f ,  r-  f )  -  H,(x,y  -  (3.64) 

The  above  equations  are  the  update  equations  for  the  x-directed  component  of  a  TE^ 
electric  field. 

4.  Derivation  of  TM^  Electric  Field  Update  Equation 

The  discrete  equivalent  of  the  TM^  magnetic  field  circulation  will  be  determined 
next.  A  sanq)le  contour  that  will  be  used  to  determine  the  electric  field  update 
equation  is  shown  below 


Hx(x,y+dy/2,t) 


Hy(x-dx/2,y,t) 


Ez(x,y,t) 


Hy(x+dx/2,y,t) 


Hx(x,y-dy/2,t) 

Figure  3.8.  A  Contour  for  TM^  Magnetic  Field  Circulation  and  Electric  Flux 

Calculations. 

The  center  of  the  surface  Sjj  is  assumed  to  have  the  coordinates  (x,  y).  We  will 
assume  a  uniform  grid  with  Ax=Ay=Al.  A  local  coordinate  system  (^,  vp)  will  be 
established,  with  the  origin  at  the  center  of  the  contour.  Any  point  (x',y’)  within  or  on  the 

contour  C^  can  be  specified  by  its  local  coordinates  ^  and  vp 
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x'=x+^  y=y  +  W  ^^^lere  -  —  <^<+—  and  -  — <v|/<+— 

The  local  coordinates  will  be  used  in  evaluation  of  the  line  and  surface  integrals  that 

constitute  the  integral  forms  of  Maxwell's  equations.  The  circulation  of  the  magnetic  field 

around  the  magnetic  field  contour  Ch  ,  assuming  counter-clockwise  reference  direction 

— ^ 

such  that  the  normal  to  the  surfece  Sg  is  in  the  direction  of  the  electric  field  E(x,y,  i),  can 
be  evaluated  approximately,  assuming  that  the  magnetic  fields  are  constant  over  the  length 

A1 ,  along  the  edges  of  the  (A1  by  Al)  square  contour 

+Hy(x'y,  07]  «  (3-65) 

[H,ix,y  -  f ,  0  +  Hyix  +  f  t)  -  H,(x,y +f,t)-  Hy(x  -  f,y,  0]a/ 

The  electric  flux  (the  "diq)lacement  current")  through  the  surface  and  the 

current  flux  (the  conduction  current)  can  be  calculated  using  the  local  coordinates 

ff  s(x'y)Ez(x'y,i)t  •dsH=  (3.66) 

JJSh 

C  O  1 T  +  ^)Ez(x + ^,y  +  V|/,  t)t  *  ~Z (£,d\y 

ff  G(x'yyEz(yy,f)z»dsH=  (3.67) 

O  G  +  +  ^)Ez{x +l,y-i-\y,()t  •~z d^d^/ 

2  2 

There  are  infinitely  many  ways  to  model  the  variation  of  the  electric  flux  denaty 

with  the  local  coordinates  ^  and  v|/  within  the  contour.  We  first  need  to  postulate  a  certain 

variation  of  the  electric  field  with  ^  and  v|/  such  that  the  integral  can  be  evaluated  over  S^. 

The  sinqjlest  model  assumes  that  the  contour  width  Al  is  small  enough  such  that  the 

electric  field  E(x+^,y+M/,t)  within  the  contour  may  be  assumed  constant  and  equal  to  the 

value  at  the  contour's  center  E(x,y,t) 

Ez(x + %,y + Ml,  0  «  Ez(x,y,  t)  (3.68) 
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This  is  equivalent  to  using  a  piece-wise  constant  approximation  of  the  actual  electric  field 
variation  with  the  z-coordinate.The  above  assun5)tion  allows  that,  although  the  electric 
field  is  constant  within  a  contour,  it  can  change  from  one  contour  to  an  other.  This  yields 
approximate  ejq)ressions  for  displacement  and  conduction  currents 

ff  <x'y)Ez(x',y',  {yz  Ezix,y,  t)  ■  <x  +  ^,y  +  \\i)d^d\y  (3.69) 

ff^  G(x'y)Ez(x'y,  i)t  •  3s  «  Ez(x,y,t)  •  0  c>(x+^,y  +  \y)i^d\\f  (3.70) 

jj  Sff  2  2 

The  integrals  of  the  permittivity  s(x-i^,y+vi/)  and  the  conductivity  a(x+^,y+vi/)  can  be 
re-written  in  the  following  manner 

^(X  +  ^,y  +  ^>)d^d^\)  =  (AO^  •  J J  s(^ + (3.71) 

T  2  L(^V  2  2 

^(x  +  ^,y  +  \\i)d^d\\i  =  (A/)^  •  J  I  jjr  a(x -i-  ^,y  +  \\f)d^d\\i  (3.72) 

T  1(^0  2  2 

The  term  in  the  brackets  is  recognized  as  the  average  permittivity  and  the 

average  conductivity  within  the  contour  Cjj. 

Savg(x,y)  =  12  0  <x+^,y+\\i)d^d\y  (3.73) 

(A/)  2  2 

Om*(x,y)  =  12  J2  ^(y+^’y  +  (3.74) 

(A/)  2  2 

The  approximate  expression  for  the  displacement  and  conduction  currents  through 
Sjj  can  now  be  written  using  the  average  permittivity  and  the  average  conductivity  as 

f[  z(x'y)Ez(x'y,  •  dsH  ~  (AySavg(x,y)Ezix,y,  t)  (3.75) 

ff  (5{x' y)Ez(x' y ,  t)t  •  dsH  « (A/)^cyovg(>!^,y)-Er(^,T»  0  (3-76) 

The  first  curl  equation  of  Maxwell  can  now  be  rqjlaced  by  an  approximate  equation 


49 


^,{[e„s(x,y)E,(x,y,  }  +  lOo^(x,y)E.(x,y,  i^KA/)’ 


wliicli  can  be  sin^lified  (because  of  media  stationarity)  to 

H:,(x,y  -  f ,  t)  +Hyix  +  f,y,  t)  -  H:,{x,y  +  f ,  t)  -  Hy^x  -  f,y,  t) »  (3.78) 

j^{E,(x,y,  /)}  •  A/  •  Savg(x,y)  +  aavg(x,y)  •  E,(x,y,  t)  •  A/ 

Tkir  equation  can  also  be  re-written  as 

§j{E,(x,y,t))  (3.79) 

— — — — AH,(x,y  -  f ,  0  +«,(x+  j,y,  /)  -  H.(y,y +",<)-  HAy  -  J.y, ')] 
AI-s^g(x,y)'- 

^avgix,y) 

The  above  equation  is  integrated  with  respect  to  time  to  get  an  approximate  equation  for 


(3.79) 


the  electric  field  at  the  present  time  t 


[^^H,<x,y-T,y)dx-^^HAx,y+^,y)<k+^^Hy{x+^,y,x)dl:-\\HAx-j,y,x)dx\ 

A  gimilar  integral  equation  can  be  written  for  the  electric  field  at  the  same  spatial 
location  but  at  an  earlier  time  t-At 

Ez{x,y,t-Lt)K  (3.81) 

A<-e^(x,joC'“  -  f.  ■')*  -  V  + f ,  T)*  +  J7"  H,<x + f  iyk- 

Subtracting  the  two  equations  we  get 
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(3.82) 


Ez(^,y,  t)  «  Ez(^,y,  t-^t)  + 


1 


AI-Sa^(x,y) 


[ll^Hz{x,y-f,T)dx-ll^Hz(x,y+f,x)dx+l_^Hyix+f,y,x)ck-j‘^_^Hy(x-f,^^^^ 


'avg\ 


The  integrals  can  be  evaluated  approximately  by  assuming  a  certain  variation  of 


magnetic  field  with  the  ten^oral  variable  t.  The  supplest  approach,  consistent  with  the 
aggiimptinns  made  for  the  field  spatial  variation,  would  be  to  assume  that  the  interval  At  is 
giTiflll  enough  such  that  the  magnetic  field  can  be  assumed  approximately  constant  within 
At  and  equal  to  the  value  at  the  center  of  the  time  interval  (t-At,  t) 


fi,(x,y  -  f.  T)  «  H,(x,y 

for 

t-At<x<t 

(3.83) 

HAx,y+f,l)oH.(x,y  +  f,l-f) 

for 

t-At<x<t 

(3.84) 

Hy(x-f,y,x)»H,(x-f,y,t-f) 

for 

t-At<x<t 

(3.85) 

for 

t- At<x<t 

(3.86) 

The  integral  of  the  electric  field  is  recognized  as  the  average  value  of  the  field 


within  the  At  interval 


Oa»g(x,y)‘At 

^avg(x,y) 


■  Ez,avg(£/)(x,y,  0 


(3.87) 


The  average  value  of  the  electric  field  is  approximately 


Ez,avg(iil)(x,y,  0  ' 


Ez(x,y,  t)  +Ez(x,y,  t-At) 


(3.88) 


The  approximate  expression  for  the  electric  field  at  location  (x,y)  and  at  the  time  t 


now  becomes 


Ez(x,y, t)  «  Ez(x,y, t-At)  + 


At 


Al-Zavg{x,y) 


(3.89) 


\HAx,y  -  f ,  /  -  f )  -  + f ,  /  -  f  >  +  f  -  f )  -  Hy(x  -  &y,  t- f )] 


gavg(x,y)  •  At  Ez(x,y,  t)  +  Ez{x,y,  t  -  At) 

^avg(x,y)  '  2 
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The  above  equation  can  be  written  in  a  somewhat  differait  form,  using  the  identity 

I  _  1  \  1  1  _  1  =  v  Z 

s  Til^VSoer  ° 


(3.90) 

vdiere  the  free  space  velocity  of  propagation  is  denoted  v^,  the  intrinsic  iicpedance  of  free 
^ace  is  denoted  Zq  (Zo  «  377Q),  and  denotes  relative  permittivity.  Introducing  the  grid 
"propagation"  velocity 

=  (3.91) 

we  can  write  the  electric  field  equation  more  compactly 

+  (3.92) 


_  7  vq  gavg(x,y)  ■  A/  E,(x,y,  t)  +  E^(x,y,  t  -  AQ 

^grid  ^r^g(p^^y)  2 

The  electric  field  "update"  equation  is  obtained  by  grouping  the  like  terms  in  the 


above  equation 


^  1  Vq  Zq  •  Oavg(x,y)  *  A/ 

~  2^^  &r,ayg(x,y) 

^  ^  1  Vo  ^0  '  ^avg(x,y)  "  ^ 


-Sr(x,y,0 


■Ez(x,y,t-At)  + 


(3.93) 


Zo 


Vo 


1 


^g’idSr^(x,y) 


.  1  Vo  Zo-Gavg(x,y)-Al 

^  ”  A  ^ 


Sr,ayg(x,y)  "  '  I'^grid  Er.mx(x,y) 

{H.(x,y-f,l-f)-  H.(x,y  +  f ,  /  -  f )  *Hy(x + f,y,  t-f)-H,(x-f,y,t-f)] 
The  equation  simplifies  for  the  case  of  non-conductive  media  (a=0) 


Vo 


Ez(x,y,  t)  «  Ez(x,y,  t- At)  +Zo— 


(3.94) 


'^g’^^^r^g(?C,y) 

[H,ix,y  -  f t  -  f )  -  H,(x,y  +  f ,  t  -  f )  +Hy(x  +  f  ,y,  t  -  f )  -  Hy(y  -  f  ,y,  t  -  f )] 

The  approximate  "update"  equations  for  TMz  magnetic  field  conq)bnents  will  be  derived 


next. 
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5.  Derivation  of  TM^  Magnetic  Field  Update  Equations 

We  will  start  with  the  Maxwell's  curl  equation  for  the  circulation  of  the  electric 


Ez{x,y, t)z  *dlE  =  z(x,y^Hz(x,y, t)x  ^Hy{x,y, t)y  ] •  •  (3.95) 

The  discrete  equivalent  of  the  magnetic  field  circulation  needs  to  be  determined 
for  two  sample  contours  a  contom  in  a  plane  parallel  to  the  yz-coordinate  plane,  and  a 
contour  in  a  plane  parafiel  to  xz-coordinate  plane. 


4 

Ez(x-Al/2,y,t^/  ^  Hy(x,y,t)^^/^ 


Ez(x  +  Al/2,y,t) 


Figure  3.9.  A  Contour  for  Updating  Hy. 

A  local  coordinate  system  (^,  Q  will  be  used,  with  the  origin  at  the  center  of  the 
contour.  Any  point  (x',2!)  within  or  on  the  contom  can  be  q)ecified  by  its  local 

coordinates 

x'=x  +  ^  z'  =  z+C  where  and 

The  circulation  of  electric  field  around  Cg  ,  assuming  counter-clockwise  reference 


direction  such  that  the  normal  to  the  surface  Sg  is  in  the  direction  of  magnetic  field 
Hy{x,y,  t),  is  therefore  given  exactly  by  the  following  sin:q)le  expression  because  there  is 
no  electric  field  variation  with  the  z-coordinate 
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(3.96) 


E,(x'y,  t)t  •  (He  = 

n E,(^x - f,y, i)z  •  d'Ct’ - E^{x  +  i)t •  di:^ 

2  2 

Eziyy,  i)z  *laE^  \Ezix  -  f,  j,  O  -  £:z(x  +  f  o]  •  a/  (3.97) 
The  magnetic  flux  through  the  surface  Sg  can  only  be  evaluated  approximately, 

since  the  exact  way  that  the  magnetic  flux  density 

By(x,y,  t)  =  \i{x,y)Hy{x,y,  t)  (3.98) 

varies  with  x  and  y  is  not  known  a  priori.  First  we  need  to  evaluate  the  flux  integral 

ILg  |t(x^/, t)Hy(x',y',  f)y  •cisE=  (3.99) 

O  n  + ^’y + + ^^y + 

2  2 

Since  the  integrand  is  not  a  function  of  the  local  z-coordinate  C  and  the  local  coordinate 

\l/=0,  the  surface  integral  simplifies  to  a  line  integral 

|r(x  +  ^,y)Hy(x  +  ^,y,t)^ •  dsE  =  A/  Q  +  ^,y)Hy(x  +  ^,y, t)d^  (3. 100) 

Se  2 

Next,  assume  that  the  contour  side  A1  is  small  enough  such  that  the  magnetic  field 
H(x+^,y,t)  within  the  contour  may  be  assumed  constant  and  equal  to  the  value  at  the 
contour's  center  H(x,y,t).  This  assumption  yields  an  approximate  expression  for  the  flux 
integral 

ff  ji,(x  +  ^,y)/fy(x+^,y,/)>^ •  dsE~ El-Hy(y,y,t)  •{_l  p(x  +  ^,y)<i^  (3.101) 

JJ  Se  2 

or,  using  permeability  averaged  in  the  x-direction 

p(x + ^,y)Hy(x  +  |,y,  t)y  •  «  (A/)'  •  Hy(x,y,  t)  ■  pa,yx,y)  (3. 102) 

JJ  Se 

The  first  curl  equation  of  Maxwell  can  now  be  replaced  by  an  approximate  equation 

[£,(X  -  f,y.  0  -  E,(x  +  f,y,  <)]  ■  a;  »  -§j{Hi^(.x.y)H,(x,y,  /))(A0"  }  (3.103) 
which  can  be  simplified  to  below  equation,  because  of  media  stationarity 


54 


Ez{x  -  f,y,  t) -Ez{x  +  ^,y,  t)  «  •  A/  •  ^{Hy(^,y,  t)  (3. 104) 


The  approxiioate  equation  can  be  re-written  as 

^  ^2,0c^J_EA.:c-f,y,O-E,(x  +  f,y,i)]  (3.105) 

The  above  equation  is  integrated  with  respect  to  time  to  get  an  approximate 
equation  for  the  electric  field  at  the  present  time  t 

~  f ■  lo 

A  similar  equation  can  be  written  for  the  electric  field  at  the  same  spatial  location 
but  at  an  earlier  time  t-At 

Hy(y,y,  r  -  AO  «  ~ 

Subtracting  the  two  equations  we  get 

Hy(x,y,  t)  -  Hy(x,y,  t-At)«  (3. 108) 


A/-|Xavg^(x,>’) 
which  can  also  be  written  as 


■ulfavilL  ^)<*  -  L  *  ^’y- 


Hyix,y,t)^  (3.109) 

Assuming  that  the  interval  At  is  small  enough  such  that  the  electric  fields  can  be 
assumed  approximately  constant  within  At  and  equal  to  the  value  at  the  center  of  the  time 


interval  At 


Ez(x  +  ^,y,x)KEz(y:  +  f,y,t-j)  for  t-At<x<t 

A/  ^  A/\ 


we  get 


Ez{x-^,y,x)KEz(^-^,y,t-i^)  for  t-At<x<t 


Hy{x,y,f)« 

H,(x,y,  t-At)  +  +  f.3'. '  -  f )  -  -  f  f)] 


(3.110) 

(3.111) 

(3.112) 
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The  term  in  the  brackets  is  recognized  as  the  average  value  of  the  electric  field 

component  at  (x,y)  within  the  interval  (t-At,  t).  The  above  equation  can  be  written  in  a 

more  compact  form  usingthe  identity 

1  _  1  1  _  _ 1 _ _  1  /  ^0  1  _  y  y  JL  ("i  I  1 

where  the  fi-ee  space  velocity  of  propagation  is  denoted  v^,  the  intrinsic  admitance  of  free 

space  is  denoted  (To  =  » 377Q)  and  denotes  relative  permeabihty. 

Zo 

Introducing  the  grid  "propagation"  velocity 

=  (3.114) 

we  can  write  the  approxumte  equation  for  the  electric  field  "update"  as 

(3.115) 

g>,>-,t-Ar)  +  y,^^^J^^_^^{£,(n4.f.;-,(-f)-£.(x-f,y.r-f)} 

In  the  case  of  free-space  (p=l),  the  equation  sinplifies  to 

Hy(x,y,t)«  (3.116) 

Hyix,y,  t-Ai)  +  Y<i^[E,ix  +  ^,y,  t-f)-E,(x-  f,y,  r  -  f )] 

The  above  equations  are  the  update  equations  for  the  y-directed  component  of  a  TM^ 
magnetic  field.  Similar  equations  will  be  formulated  next  for  the  x-dire(rted  magnetic  field 
component. 
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A 


^/Ezj(x,y+dl/2,t) 


— I  Hx(x,y,t) 


Ez(x,y-dl/2,t) 


Figure  3.10.  A  Contour  for  Updating  E[^. 


A  local  coordinate  system  will  be  used,  with  the  origin  at  the  center  of  the 
contom.  Any  point  (y’jZ')  within  or  on  the  contour  can  be  specified  by  its  local 


coordinates 


y' +  \^diere  — ^<v|/<+^  and  — 

By  following  the  same  procedure  shown  for  update  equation,  we  can  determine 

the  approximate  equation  for  the  magnetic  field  "updates'*  as 

H,(x,y,t)>^  (3.117) 

Hx(x,y,  t-M)  +  ^0  "  f’  ^ "  f ) +  f’  ^  "  f )  } 

In  the  case  of  fiee-space  (|t=l),  the  equation  sin:q)lifies  to 

H.(x,y,t)«H,(x,y,t-At)  +  Yo^[E,(x,y-^,t-f)-E^x,y*^,t-!^)\  (3.118) 

The  above  equations  are  the  update  equations  for  the  x-directed  conq)onent  of  a  TE^ 
electric  field. 
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B.  TRANSPARENT  GRID  TERMINATION  IN  2-D 


1.  2-D  TE^  and  TM^  Grids 

The  electric  and  magnetic  field  update  equations  have  been  derived  using  local 
coordinate  systems,  with  origins  at  the  centers  of  the  magnetic  and  electric  contours, 
respectively.  These  equations  now  need  to  be  "converted"  to  a  global  coordinate  system, 
that  is  to  the  grid  of  equi-distant  saropling  points  in  the  xy-plane.  Since  we  have  either 
TE^  or  TM^  2-D  electromagnetic  fields  there  will  be  TE^  and  TM^  grids  as  well.  Let  us 
assume  that  our  domain  is  an  L  by  L  square.  We  will  assume  a  uniform  grid,  that  is  a 
discretization  step  A1  that  is  not  changing  with  position.  The  fields  are  then  "sattq)led" 
nging  a  spatial  step  A1=L/(N-1).  The  electric  and  magnetic  field  san^le  locations  are 
"interleaved",  as  shown  below  for  N=5.  The  spatial  grid  edge  saioples  can  be,  in  princ^le, 
either  electric  field  sauries  or  magnetic  field  samples.  We  will  select  electric  fields  for 
spatial  edge  sanoples  for  TE^  grids  and  magnetic  fields  for  TM^  grids.  The  "generic"  2-D 
grid  below  is  thus  applicable  to  either  TE^  or  TM^  fields. 
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Figure  3.11.  2-D  Grid 

An  (N  by  N)  TE^  grid  has  (NXN)  magnetic  field  nodes,  (N-2XN-1)  electric 
field  nodes,  and  (N-l)(N-2)  E^  electric  field  nodes.  Similarly,  an  (N  by  N)  TM^  grid  has 
(N-2)(N-1)  E^  electric  field  nodes,  (N-lXN-2)  magnetic  field  nodes,  and  (N+lXN) 
magnetic  field  nodes.  TE^  electric  and  magnetic  field  "update"  equations  can  now  be 
"converted"  to  electric  and  magnetic  field  grid  "update"  equations  by  replacing  the 
variables  x,  y  and  t  with  the  grid  coordinates  of  the  electric  field  ^atial  and  tenq)oral 
sampling  points 

X  iM,  i  =  0,  L..N  y  JM,  i  =  0,  \...N  t  ->  nAr,  «  =  0, 1..JV, 
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The  notation  can  be  sunplified  by  omitting  the  common  A1  and  At  terms,  and  using 


a  superscript  for  the  index  of  the  ten^oral  sanq)ling  point.  The  grid  equations  for 


conductive  media  are  obtained  from  equations  (3.62 , 3.65) 

J  1  Vo  ^0  •  A/  •  o avgy(i,j) 

EV.UJ)  - 

J  _j_  1  *  ^avgy\}:,J) 

2  S  r^avgy 


(3.119) 


-  - — j--  Hz  (i,j  +  -^)-Hz  (/,y-2) 

j  _j_  1  Vq  '  ^avgyK^J)  L 

2  ^ grid  £  r,avgy(J’  <ij) 

J  1  ^0  ^0  *  '  ^avgxi}->j) 

^r,avgx(i,j)  pr>-\(i  A 

J  ^  1  Vq  Zq  •  a/  •  OgygxQj) 

2  Sr,avg^(j^j) 

rx  Vo  1 


(3.120) 


1  Vo  ^0  '  *  ^avgx\hj)  L 


- ^ 

2  ^r^avgxif^jf) 


Vq  1 


(3.121) 


Ex'ku - i) -£r"(/,y +  i  t)  +Ef\i  +  lj)-EfHi-^J) 


The  TE  electric  and  magnetic  field  grid  update  equations  for  non-conductive 


media  (cj=0)are 


E^(iJ)  « ^ O' J) + ^0 ^(F^  ‘  +  '(V-l)  (3.122) 

^;(/,y)  «£rH^y)  +  <3.123) 


(3.123) 
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f;(v)  =  HT'  (ij)  +  foTST— ^ 


^gnd  \ir,avg(j,j) 


(3.124) 


eTHu + i)  -  eTHu- i) -  hJ) -Ef^  0 + \J) 


Finally,  the  free-space  TE^  grid  equations  are 


E"MJ)  -  E;-'  (i,J)  +  z„^ .  -  \J)  -  Hf\i  + 1  j) 


(3.125) 


(3.126) 


(3.127) 


HT\iJ)  +  Yo^  •  Ef\iJ+\)-Ef\iJ-\)+Ef\i-\j)-E7\i^\,j) 


Dual  TM^  grid  equations  for  conducting  media  are 


+  (3.128) 

(3.129) 


j  1  Vq  Eq  •  0’avg(/,y)  •  A/ 
2  ^r,avg(j,J^ 

1  Vq  Zo-Oavg(iJ)  Al 
2^gnd  ^r,avg(i,j) 


ET-\hJ)+ 


(3.130) 


Vo  1 


J  ^  1  ^0  ^0  *  ^avg\hj)  '  L  -» 

2  Zr,ctvgi}^f) 

The  TM^  electric  field  grid  equation  sinq)lifies  for  non-conductive  media 


E”iiJ)«E”  \i,j)+ZQ 


Vq  1 
^S'^^Zr^gQj) 


(3.131) 


Finally,  the  free-space  TMz  grid  equations  are 

ffiO'J)  =  HT'  (V)  +  I'o^^  ■  \eT'\i,J  -\)-  Ef‘  (iJ + j) 


(3.132) 
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(3.133) 


E'^,iiJ)«ErKiJ)+Z,^^- . 

Hl'kiJ  -  \) -HTkiJ  +  i)  +Hf\i  +  \J)  -  h";'  (/•  -  \J) 


(3.134) 


Note  that  TE^  electric  and  magnetic  field  grid  equations  in  2-D  have  the  following 
general  form 

Er  =  CExxEf  +  CE^^^Hf  (3.135) 

-  =  CEyxEf  +  CEya'^Hf  (3.136) 

Hr  =  CmHf +  Cim^Ef +  CHyi^Ef  (3.137) 

where  C's  are  constants  (real  numbers)  that  depend  on  the  media  properties  and  the 

velocity  ratio  \Jv  ^  ,  and  "del"  operator  represents  the  spatial  derivative  (gradient). 


Similarly,  the  TM^  electric  and  magnetic  field  equations  can  be  written  as 

Hr  =  CH.iHf +CH:a'^Er 

Hr  =  CnyiHf  +  CHyi^Ef 
Er  =  CExEf  +  CE:a'^Hr  +  CEyi'^Hf 


(3.138) 

(3.139) 

(3.140) 


2.  2-D  Grid  Termination 


The  grid  equations  derived  previously  are  valid  for  all  the  nodes  except  for  the 
nodes  on  the  grid  edges.  The  reason  is  that  a  grid  edge  node  for  a  field  component  in  a 
transverse-to-z  plane  has  only  one  neighbor  node  instead  of  four  like  a  node  that  is  not  on 
the  grid  edge,  as  shown  below. 
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> 


"made"  node  s 


"just  inade"  boundary 


"edge"  node 


boundary 


Figure  3.12.  Edge  Nodes  of  a  2-D  Grid. 

An  equation  for  the  circulation  of  the  z-directed  field  conqjonent  can  not  be 
written  for  an  edge  node  because  a  z-directed  field  con:q)onent  on  the  edge  shown  in  the 
figure  above  can  not  be  updated  using  the  grid  equations  discussed  so  fer.  Extending  the 
ideas  of  transparent  grid  termination  (TGT)  and  the  discrete  boundary  impulse  response 
(DBIR)  applied  for  1-D  problems,  the  update  equations  for  edge  nodes  can  be  devised 
based  upon  convolutions  of  the  fields  one  layer  inward  and  the  pre-calculated  impulse 
responses  from  the  inward  layer  to  the  grid  edge  nodes.  An  edge  field  will  be  expressed  as 
a  superposition  of  responses,  as  illustrated  in  the  figure  below  for  a  top  edge  grid  node. 
(The  same  apphes  to  other  edge  nodes  as  weDL) 
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Figure  3.13.  Edge  Field  as  a  Superposition  of  Responses. 

The  superposition  can  be  also  "visualized"  using  the  concept  of  a  multiport .  The 
nodes  on  grid  edges  can  be  considered  as  output  ports  of  a  linear  multi-port.  Assuming, 
for  simplicify  that  the  grid  is  square  with  N  edge  field  nodes  on  each  ride,  the  total 
number  of  output  ports  will  be  4(N-1).  The  nodes  on  the  layer  just  inside  the  grid 
boundary  will  be  considered  as  input  ports  of  the  linear  multiport.  There  will  be  a  total  of 
4(N-3)  input  ports.  The  equivalent  multiport  is  shown  below. 


Figure  3.14.  Modeling  2-D  Grid  Boundary  as  a  Multiport. 
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The  input  and  output  fields  can  be  either  electric  or  magnetic  fields.  Furthermore, 
the  input  and  output  fields  need  not  be  of  the  same  "kind",  that  is  the  input  fields  could  be, 
in  principle.  E-fields  and  the  output  fields  could  be  H-fields,  or  vice  versa.  However,  we 
will  select  the  input  and  output  fields  to  be  of  the  same  "land",  that  is  either  all  £  or  all  H 
fields,  such  that  the  in^ulse  responses  h^(t)  will  be  "dimensionless".  The  input  and 
output  fields  will  be  H-fields  for  TE^  and  E-fields  for  TM^  fields.  In  this  manner,  the  same 
impulse  responses  h^(t)  can  be  used  for  both  TE^  and  TM^  grids.  The  equations  will  be 
derived  for  TE^  fields,  that  is  with  electric  fields  at  the  grid  edges.  Dual  equations  for  the 
TM^  fields  can  be  obtained  surply  by  replacing  E's  with  Hs  in  the  TE^  equations. 

There  are  four  sides  of  the  square  grid  boundary.  The  top  side  has  the  electric 
field  E^(0,y,t),  the  bottom  side  has  the  electric  field  E^(L,y,t),  the  left  side  has  the  electric 
field  E/x,0,t),  and  the  right  side  has  the  electric  field  E/x,L,t).  The  layer  just  inside  will 

have  the  corresponding  fields  as  E^(Al,y,t),  E^(L-Al,y,t),  E^(x,Al,t),  and  E^(x,L-Al,t).  The 

limits  on  the  values  of  x  and  y  for  the  edges  are  0  and  L  (0  <  x,y  <  L),  while  the  limits  fi>r 
the  layer  just  inside  are  A1  and  L-Al  (A/  <  x,y  <L-  A/).  Ushig  the  grid  coordintes  (ij)  the 

fields  on  the  grid  boimdary  will  be,  starting  at  the  upper  left  comer  (coordinates  0,0)  and 

using  a  clockwise  reference  directionE"(0,y), £’"(/,  l,y),£"(/,0).  These  can 

be  grouped  into  a  "vector"  of  electric  fields  on  the  grid  boundary 

E,ounaary(i)=^{EMj,t)  E,{N-\J,t)  E,{i,0,t)f  (3.141) 

where  T  indicates  transposition  (the  vector  is  a  column  vector,  but  it  has  been  written  as  a 

transpose  of  a  row  vector  to  save  space).  In  an  analogous  manner  an  electric  field  vector 

can  be  formed  for  the  layer  just  inside  the  grid  boundary 

Ejustjnsiae(i)  =  [E,{lJ,t)  E,(i,N-2,t)  E,{N-lJ,t)  E,{i,\,f)f  (3.142) 
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The  advantage  of  the  vector  arrangement  is  that  we  can  now  use  a  single  index 
(say  index  m)  to  identify  a  field  on  the  grid  boundary  and  another  single  index  (say  index 
n)  to  identify  an  electric  field  on  the  layer  just  inside  the  grid  boundaiy.  The  use  of  single 
indices  m  and  n  allows  for  a  direct  and  contact  refefrence  to  the  multiport  rq)resentation 
of  the  grid  termination  problem  The  electric  field  on  the  grid  boundary  Eb„^j^(t)  can  be 
expressed  as  a  convolution  of  the  electric  fields  at  the  layer  just  inside  the  grid  boundary 

^just  insideC^)  *  matrix  of  iir^ulse  responses  H(t). 

E bounda?y(t')  ”  -^(0  ^  (3.143) 

where  the  4(N- 1)  by  4(N-3)  matrix  of  in^ulse  responses  can  be  written  as 

hia(t)  /2u(0 . *i,4(jv-3)(0 

h2,l(i) . h2,4(N-3)(t) 

H(t)=  .  (3  144) 


[  h4(N-l)At)  hA(N-\)Ai) . ^4(Ar-l),4(V-3)(0  J 

The  electric  field  at  a  boundary  node  (identified  by  a  particular  value  of  the  index 

m)  can  be  written  as  the  product  of  the  m-th  row  of  the  in^ulse  response  matrix  H(t)  and 

the  column  vector  of  the  electric  fields  just  inside  the  boundary  (indentified  by  index  n) 

Eboundaryim,  0  ~  (3"  145) 

^^V-3)  ^  Ejus,_inside(n,  t)  =  £  Ejustjmidein,  t)  •  t)^* 

v^iiere  the  summation  is  over  all  the  nodes  on  the  layer  just  inside  the  boundary.  The 
discretized  forms  of  the  above  equations,  involving  san:q)les  of  electric  fields  at  t=kAt,  are 
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the  gims  of  discrete  convolutions.  Since  a  discrete  convolution  is  itself  a  sum,  the  result 
for  the  boundary  nodes  will  involve  double  summations  (the  expression  below  uses 
superscripts  for  time  san^les  and  convolution  summation). 

(«)  ■  <3-  W6) 

The  convolution  equations  express  the  fields  at  the  edges  as  weighted  sums  (the 
weighting  coefficents  are  the  values  of  the  discrete  boimdary  impulse  responses  h^)  of  the 
time  histories  of  the  fields  "just  inside"  the  grid.  In  that  respect  they  are  equations  of  the 
same  type  as  the  equations  for  the  non-edge  grid  nodes  (equations  ??  and  ??),  except  that 
they  involve,  in  general,  summations  with  many  more  terms.  However,  the  impulse 
responses  hj„^(t)  converge  to  zero  relatively  fest  which  effectively  reduces  the  munber  of 
significant  terms  in  the  convolution  sums.  The  convolution  sums  can  thus  be  truncated 
such  that  only  a  certain  munber  of  the  most  recent  values  of  the  electric  fields  just  inside 
the  boundary  are  needed  (this  reduces  the  number  of  electric  fields  that  have  to  be  stored 
and  updated).  The  discrete  boundary  in:q)ulse  responses  h^(kAt)  need  to  be  determined 
(and  saved)  only  once,  just  like  in  1-D,  for  a  selected  grid  velocity  Vgjij=Al/At.  The 
impulse  reqronses  can  be  stored  in  a  rectangular  matrix  form.  A  particularly  suitable 
matrix  form  would  have  4(N-1)  rows,  one  row  per  boundary  node,  and  4(N-3)N^ 
colunms,  wiiere  denotes  the  "truncated"  length  of  the  discrete  boimdary  impulse 
responses.  This  form  is  suitable  because  the  boundary  electric  fields  can  then  be 
calculated  as  a  product  of  such  a  matrix  and  a  column  vector  constracted  of  4(N-3) 
subvectors  of  length  N^.  The  subvectors  represent  the  most  recent  values  of  the 
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electric  field  at  nodes  just  inside  the  boundary.  (The  choice  of  N,j  depends  on  the  desired 
arqphtude  "resolution",  with  the  values  over  20  giving  very  good  results  for  Vgrij/vo  =  JD 
,  D=2  for  two  dimensional  problems).  These  subvectors  are  updated  at  each  time  step  by 

shifting  all  of  their  values  down  (that  is  "fiither  into  the  past")  by  one,  discarding  the  last 

value,  and  updating  the  first  subvector  element  with  the  most  recent  field  value  calculated 

via  standard  grid  update  equations.  The  double  summations  are  thus  efficiently  executed 

as  matrix-vector  multipHcations. 

The  discretized  in:q)ulse  responses  h-^  are  determined  in  a  maimer  analogous  to 
the  1-D  case,  using  the  discrete  equivalent  of  the  delta  fimction  which  we  will  denote  as 
5'‘.  It  is  convinent  to  use  the  equivalent  multiport  to  explain  the  procedure.  In  order  to 
determine  the  discrete  boundary  impulse  response  for  a  particular  input  port  n,  all  input 
ports  except  the  n-th  input  port  are  set  to  zero  and  delta  impulse  is  appHed  to  the  n-th 
input  port.  Recording  aU  the  outputs  fi-om  t=0  to  t=NtAt  provides  us  with  4N  discrete 
boundary  impulses  of  duration  N^. 

=  (3M7) 

=  ^  1...4(Af-3)  and  p*n  (3.148) 

The  process  is  repeated  4(N-3)  times  (since  there  are  4(N-3)  input  ports)  and  the 

results  are  stored  in  a  4(N-1)  by  4(N-3)Nb  matrix.  This  matrix  is  then  used  to  inqilement 
the  transparent  grid  termination  via  the  matrix  multiplication  by  the  4(N-3)N^  column 
vector  of  the  time-histories  of  the  fields  just  inside  the  boundary.  Depending  on  the  values 
of  Nh  and  N,  there  will  be  input  nodes  that  are  sufficiently  fer  from  an  output  node  such 
that  the  time  it  takes  for  an  input  impulse  to  propagate  to  the  ouput  port  exceeds  Nj,.  In 
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such  a  case,  the  contribution  of  such  an  input  node  to  the  output  node  is  known  to  be  zero 
in  advance.  This  may  be  used  to  reduce  the  number  of  operations  in  a  TGT 
implementation.  The  process  of  obtaining  the  2-D  discrete  boundary  impulse  responses  is 
carried  out  on  a  "large"  grid  whose  size  should  be  at  least  (N+N^  by  N+NJ,  or 
equivalently,  the  distance  from  the  TGT  boundary  to  the  boundary  of  the  "large"  grid 
should  be  at  least  N^Al/2.  The  termination  of  this  "larger"  grid  (in  typical  applications 
N»N^  and  the  "larger"  grid  would  be  only  incrementaly  larger)  is  immaterial,  since  the 
reflections  off  its  boundary  do  not  arrive  at  the  output  nodes  (where  they  would  have 
corrupted  the  impulse  responses)  before  the  time-stepping  has  been  terminated.  When 
calculating  the  discrete  boundary  iirqrulse  responses  (DBIR),  the  nodes  interior  to  the 
TGT  bormdary  need  not  be  updated  since  all  the  nodes  on  the  layer  just  inside  the  TGT 
boundary  (the  input  ports)  are  zero  for  t>0  (at  t=0  only  one  node  on  the  layer  just  inside 
the  TGT  boundary  has  the  value  of  1).  The  DBIR  calculations  thus  require  updates  only 
for  the  nodes  between  the  TGT  boundary  and  the  large  grid  boundary  which  can  result  in 
significant  computational  time  savings  if  N»N^.  The  update  equations  used  for  the 
region  between  the  "iimer"  (TGT)  and  the  outer  boundaries  are  the  free  q)ace  equations, 
smce  this  region  is  assumed  to  be  free  space.  (The  medium  between  the  two  boundaries 
can  be  any  homogenous  medium  extending  to  infinity,  but  the  most  practical  case  is  the 
free  space.) 
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Figure  3.15.  2-D  Grid  for  Calculation  of  TGT  Boundary  Inqjulse  Responses. 

The  shape  of  the  discrete  boundary  impulse  re^onse  h^(t)  determined  as 
described  in  this  section  will  depend  on  the  relative  positions  of  the  observation  (output, 
index  m)  and  the  source  (input,  index  n)  nodes.  In  general,  the  more  the  indices  diflfer 

from  each  other  the  following  will  be  noticed 

•  the  impulse  response  will  start  later; 

•  the  maximum  value  of  the  impulse  response  will  be  smaller; 

•  the  impulse  response  will  be  more  "spread  out"  in  time. 

To  illustrate  these,  three  sample  impulse  responses  are  shown.  The  observation 
point  (output  port)  is  the  same  for  all  three  responses,  and  is  positioned  in  the  center  of 
one  side  of  the  TGT  boundary.  The  three  discrete  botmdaiy  impulse  responses  are  the 
responses  to  the  delta  function  applied  to 


•  a  source  (input  node)  immediatelly  below  the  output  node; 

•  an  input  node  5  nodes  to  the  side  of  the  source  node  immediately  below  the 
output  node; 

•  an  input  node  10  nodes  to  the  side  of  the  source  node  immediately  below  the 
output  node. 

The  velocity  ratio  for  the  DBIR  calculations  was  set  to  l/Jl.  The  plotting 
program  interpolates  linearly  between  the  saropling  points  kAt.  The  impulse  response 

were  truncated  to  a  duration  of  Nh=40.  Because  of  the  vast  change  in  scale  (the  largest 

value  of  DBIR  for  a  source  node  just  below  is  1/2  v^Me  the  largest  DBIR  value  for  a 

sotuce  5  nodes  to  the  side  is  about  25  tunes  smaller,  and  the  largest  DBIR  value  for  a 

source  node  10  nodes  to  the  side  is  about  1,000  times  smaller)  the  in5)ulse  req)onses  are 

shown  mdividually.  The  rapid  decrease  of  the  impulse  req)onse  peak  values  with  distance 

offers  the  possibility  to  implement  TGT  using  "local"  equations  instead  of  using  "global" 

convolutions.  An  approximate  "local"  TGT  inq)lementation  would  involve  only  a  few 

input  nodes  nearest  to  an  output  node,  \^Me  a  "global"  TGT  in:q)lementation  mvolves  all 

input  nodes  (aU  the  nodes  on  the  layer  just  inside  TGT).  The  loss  of  accuracy  due  to 

"localization"  is  not  excessive,  because  of  the  rapid  reduction  of  the  maxima  of  h„,^(t)  as 

the  difference  between  m  and  n  iacreases.  Fnially,  the  DBIR  shown  are  equally  applicable 

to  TE^and  TM^  grids. 
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Figure  3.17b.  DBIR  For  A  Source  Node  10  Nodes  to  the  Side 

C.  2-D  TGT  RESULTS 

Li  order  to  test  the  TGT  we  will  use  the  TMj  grid  model  shown  in  Figure  3.11. 
Therefore,  the  nodes  on  the  grid  edges  are  E-field  nodes.  Note  that  the  results  would  be 
the  same  for  TE^  fields,  with  H-fields  node  on  the  edges  of  the  grid.  The  source  is  a 
unit-anq)litude  delta  inq)ulse  of  E^  electric  field  applied  at  the  center  node.  The  inq)ulse 
source  at  the  center  of  the  grid  gives  rise  to  a  cylindiical  wave  that  propagates  radialfy 
away  from  the  grid's  center.  The  cylindrical  wave  is  not  perfect,  as  its  is  "distorted"  by  the 
grid  and  develops  a  "wake"  as  it  propagates  through  file  grid.  The  cylindrical  wave  is 
incident  to  the  TGT  boundary.  If  the  TGT  were  pofect,  the  total  power  within  the  grid 
would  be  the  same  as  the  total  power  within  the  same  region  for  an  infinite  grid. 
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However,  since  the  TGT  is  not  perfect  (because  of  truncaton  of  impulse  responses)  there 
will  be  some  increase  in  the  total  power  that  comes  from  the  "reflection"  off  TGT.  A 
measure  of  TGT  performance  in  2-D  will  be  the  relative  increase  in  the  total  power,  that  is 
the  ratio  of  the  power  increase  for  the  TGT  and  the  power  for  an  infinite  grid.  This  ratio 

is  a  function  of  time  and,  since  it  is  normally  very  small  it  is  best  represented  in  dB 

^  1  -  (Pm[-PTGT\ 

Qtgt=  10-logjo|^ — J 

We  have  denoted  this  ratio  as  since  it  e)q)resses  the  TGT  "quality".  We  will 
select  a  square  grid  with  21  nodes  on  each  side  for  our  2-D  TGT  test.  The  inpulse 
responses  will  be  truncated  such  that  the  TGT  performance  can  be  observed  for  various 
selected  impulse  response  durations.  Note  that  the  power  within  the  grid  will  be  the  same 
for  the  case  of  an  infinite  grid  and  the  TGT  until  the  outward-propagating  cylindrical  wave 
reaches  the  TGT.  All  results  will  be  shown  for  the  case 

"^grid  —  ■  Vo 

The  grid  power  and  the  TGT  quahty  fector  0^^  will  be  plotted  as  functions  of 
time,  with  the  coordinate  origin  (t=0)  reffering  to  the  time  when  the  outgoing  cylindrical 
wave  first  reaches  the  TGT. 

First,  we  truncate  the  impulse  responses  such  that  their  duration  is  t^=  lOAt.  The 
power  within  the  grid  is  plotted  as  well  as  the  power  difference  and  the  relative  power 
difference  in  dB  (the  Q^ct).  The  plot  of  power  vs  time  ^ows  the  decrease  in  power 
towards  0  as  the  cylindrical  wave  leaves  the  grid.  If  the  grid  and  the  TGT  were  ideal,  the 
power  will  be  constant  until  the  wave  reaches  the  grid  nodes  at  the  centers  of  each  side 
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(the  outgoing  wave  reaches  these  nodes  first  since  they  are  closest  to  the  grid  center 
wliere  the  wave  has  originated).  The  power  will  then  decrease  towards  zero  as  more  and 
more  of  the  wave  leaves  the  grid.  Finally,  the  wave  would  reach  the  four  comers  of  the 
grid  (these  nodes  are  the  farthest  firom  the  grid  center  and  the  wave  gets  to  them  last)  and 
leave  the  grid  alltogether,  the  power  within  the  grid  being  zero  firom  then  on.  However, 
since  neither  the  grid  nor  the  TGT  are  perfect,  there  will  he  some  small  "residual"  energy 
within  the  grid  at  all  times,  converging  to  zero  as  time  progesses.  The  relative  power 
difference  in  db  (the  Q^ct)  curve  shows  that,  for  t<  t^  the  power  "reflected"  off  the  TGT  is 
about  150  dB  below  the  power  that  the  same  domain  (within  the  TGT  boundary)  would 
have  for  an  infinite  grid.  For  t  >  t^  the  power  "reflected"  off  the  TGT  is  more  that  20  dB 
below  the  power  of  an  infinite  grid.  This  means  that  the  effects  of  the  TGT  are  about  two 
orders  of  magnitude  smaller  than  those  of  the  grid  "imperfection"  for  t  >  tj,  and  about  15 
orders  of  magnitude  smaller  for  t  <  t,,.  TTie  plots  for  t^=20, 1^=30,  and  1^=40  confirm  the 
above  as  welL  Note  that  increasing  the  duration  of  the  impulse  response  also  reduces  the 
"reflections"  off  TGT  for  t>th.  Also,  note  that  for  a  grid  of  N  nodes  on  the  ade  the 
duration  of  in^ulse  responses  greater  than  2N  allows  for  any  node  on  the  layer  just  inside 
TGT  to  contribute  at  least  one  value  to  any  of  the  nodes  on  the  TGT  (because  there  are 
2N  nodes  fi:om  one  comer  of  the  TGT  to  the  opposite  comer). 
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Figure  3.18.  Residual  Power  for  =3At. 


Residual  Power  for  TGT  Boundary 
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Residual  Power  Difference  forTGT  &  infinite  grid 
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Figure  3.19.  Residual  Power  Difference  Between  TGT  and  Infinte  Boundary  for  t(,-3At. 
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Figure  3.20.  Residual  Power  for  t^=10At. 


79 


Residual  Power  Residual  Power 


Sample  Times  Since  Touching  TGT  grid 


Figure  3,22«  Residual  Power  for  t^— 20At. 


81 


82 


Linear  Scale 


Log  Scale  (dB) 


Duration  of  Impulse  Response 


Figure  3.28.  Maximum  Power  Difference  Between  TGT  and  Infinte  Boundary  as  a 
function  of  lirqjulse  Response  Duration  (t^). 
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Residual  Power  for  TGT  Boundary  for  vr=1/[2*sqrt(2)] 


Sample  Times  Since  Touching  TGT  grid 


Figure  3.30.  Readual  Power  for  t^=40At,  using  Vr  = 
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Power  Difference: Log  Scale 


Sample  Times  Since  Touching  TGT  grid 


F^itre  3o31.  Residual  Power  Difference  Between  TGT  and  Infinte  Boundary  for  t,,-40At, 

using 
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IV.  ANALYSIS  OF  TGT  FOR  3-D  FD-TD 


A.  FD-TD  FORMULATION  DJ  3-D 

1.  FD-TD  Equations  in  3-D 

The  incident  and  the  scattered  electromagnetic  fields  and  the  media  parameters  in 
3-D  problems  vary  with  three  spatial  coordinates(x,y,z)  ,  so  the  electric  and  magnetic 
fields  will  be  fimctions  of  the  spatial  coordinates  x,  y,  z,  and  time  t.  A  3-D  electric  field 
vector  E  can  be  represented  in  terms  of  its  othogonal  components  as  shown  below 


Ex 


. - . >  y 


Figure  4.1.  Electric  Field  Vector  Components. 

'E(x,y,z, t)  =  Ex(x,y,z, t)x  +Ey(x,y,z, t)y  +Ez(x,y,z, t)  z  (4.1) 

The  direction  of  propagation  is  indicated  by  the  direction  of  the  velocity  of 

propagation  vector  ~v  .  Similarly,  for  the  magnetic  field  vector  H  we  can  write 

~H(x,y,z, t)  =  H:,ix,y,z, t)x  +Hy(x,y,z, t)y  +Hz{x,y,z, t)  z  (4.2) 

The  electric  and  magnetic  fields  are  solutions  of  Maxwell's  equations 
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(4.3) 


(4.4) 


dt 


Jls;,  •  +1/5^  a(x,y,z)E(x,y,z,  t)  •  dsn 


The  second  curl  equation  does  not  have  the  source  current  term,  since  it  has  been 
assumed  that  there  were  no  source  currents  in  the  domak  of  interest.  The  contours  of 
integration  for  the  electric  and  the  magnetic  field  circulations  are,  in  general,  different  and 
will  thus  labeled  for  an  electric  field  contour  and  Cjj  for  a  magnetic  field  contour. 
Similarly  the  surfaces  associated  with  the  contours  will  be  labeled  Sg  for  the  magnetic  flux 
and  S„  for  the  electric  flux.  In  order  to  determine  the  discretized  forms  of  Maxwell's  curl 
equations  in  the  integral  form  for  the  3-D,  we  will  use  3-D  Yee  electric  and  magnetic  cells 
[Ref  1].  One  such  cell  is  shown  below. 
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The  electric  field  contoixrs  Cg  and  the  magnetic  field  contours  are  square  (A1  by 
Al)  contours.  The  linked  Cg  and  contours  are  aU  orthogonal  to  each  other,  evoking  the 
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idea  of  electric  and  magnetic  field  linkage  in  all  directions.  Based  on  Yee's  cells,  the 
electric  and  magnetic  contours  and  their  associated  surfaces  will  be  used  to  discretize 
Maxwell's  curl  equations. 

2.  Derivation  of  Magnetic  Field  Update  Equations 

We  need  to  derive  the  update  equations  for  each  of  the  three  components  of 
magnetic  field  vector  H.  A  contour  for  deriving  the  update  equation  for  H,j(x,y,z,t)  is 
drown  below.  This  contour  can  be  thought  of  as  the  firont  face  of  the  Yee  cell  (cube)  from 
the  previous  figure.  Note  that  the  coordinates  of  the  four  electric  field  components  are 
expressed  in  terms  of  the  coordinates  of  the  contour’s  center. 


The  local  coordinates  will  be  used  in  the  evaluation  of  the  line  and  surface  integrals 

that  constitute  the  integral  forms  of  Maxwell's  equations.  The  circulation  of  the  electric 

field  around  the  electric  field  contour  ,  assuming  counter-clockwise  reference  direction 

such  that  the  normal  to  the  surface  is  in  the  direction  of  the  magnetic  field  Hx(^,y,z,  t), 
can  be  evaluated  approximately,  assuming  that  the  electric  fields  are  constant  over  their 

re^ective  countoufs  sides  of  length  A1 

Ey(x'y,z', t)x  +E^(x',y',z', t)y  ]  •  «  (4.5) 

Ey(x,y, z-^,t)-  Ey(x,y,  ^  +  y,  0  +Ez(x,y  +  z,  f)  -  E,(x,y  -  z,  t)  •  A/ 

The  magnetic  flux  through  the  surface  (the  right-hand  side  of  the  first  curl  equation) 
can  be  calculated  using  the  local  coordinates 

li(x',y', z')Hz(x'y, z', f)z  •dsE=  (4.6) 

JJ  JJ  \i.{x,y  +  \\i,z  +  (^H:c{x,y-\-\^,z  +  <;,t)x  •  xd\\idq 

2  2 

There  are  infinitely  many  ways  to  model  the  variation  of  the  magnetic  flux  density 

(with  reject  to  local  coordinates  v|;  and  Q  within  the  contour.  We  first  need  to  postulate 

a  certain  variation  of  the  magnetic  field  with  v|/  and  ^  such  that  the  integral  can  be 

evaluated  over  S^.  The  simplest  model  assumes  that  the  contour  width  A1  is  small  enough 

such  that  the  magnetic  field  H(x,y+v|/,z+^,t)  within  the  contour  may  be  assumed  constant 

and  equal  to  the  value  at  the  contour’s  center  H(x,y,z,t). 

Hx(x,y  +  y\i,z+q,t)!>iH:c(x,y,z,t)  (4.7) 

This  is  equivalent  to  using  a  piece-wise  constant  (or  2-D  "pulse"  expansion) 

approximation  of  the  actual  magnetic  field  variation  with  respect  to  y  and  z.  Note  that  the 
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above  assun^tion  allows  that,  although  the  magnetic  field  is  constant  within  a  contour,  it 
can  change  from  one  contour  to  an  other.  This  yields  an  approximate  expression  for  the 
for  the  magnetic  flux 

ff  p(x'y,z')H,ix',y',z',t)z  •  'ZE«H,{x,y,z,t)-p^  f J  \x{x,y+yy,z  +  q)d\ydq  (4.8) 

JJ  5^  2  2 

The  integral  of  the  permeability  |i(xy+\|/,z+C)  can  be  re-written  in  the  following  manner 

r!  rt  + M/,  2 = (A/)^  •  n  n  ^ 

J  J  __  /  A  /\  V  * 


,y  +  \]i,z  +  yi\\id(;\  (4.9) 


The  term  in  the  brackets  is  recognized  as  the  permeability  averaged  in  the  y 
and  z-directions  within  the  contour  Cg  (vvfiich  is  in  a  x  =  const  plane) 

Via.gxix,y,  z)  =  0  +  ^,z  +  <;)dwd<;  (4. 10) 

(A/)  2  2 

The  approximate  expression  for  the  magnetic  flux  through  Sg  can  now  be  written  using  the 
average  permeability  as 

li(x'y,z')H:,(x'y,z',t)x  •dsE~  {Mf\Jia^gx{x,y,z)Hxix,y,z,t)  (4.11) 
Again,  this  approximate  expression  resulted  from  the  piece-wise  constant 
approximation  of  the  magnetic  field  with  respect  to  y  and  z-coordinates.  The  main 
advantage  of  the  piece-wise  constant  expansion  enq)loyed  above  is  its  afr5)licity.  The  first 
curl  equation  of  Maxwell  can  now  be  replaced  by  an  approximate  equation 

Ey(x,y,  2  -  y ,  0  -  Ey{x,y, ^  +  yO  +  E,{x,y  +  y , z,  0  -  Ex(x,y  -  z,  t)  «  (4. 12) 

~{[\ya^gx{x,y,z)Hx{x,y,z,  0]A/} 

which  can  be  simplified  (because  of  media  stationarity)  to 

Ey(x,y,z-^,t)-Ey(x,y,z^^,t)+E,(x,y^^,z,{)-Ex(x,y-^,z,f)»  (4-13) 

-^{H^{x,y,z,{)}  ■  M •  \Xavgx{x,y,z) 
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This  equation  can  ako  be  re-written  as 


(4.14) 


77 - 7 - r  \Eyix,y,z-^,()- Ey{x,y,z+  /)  +E,{x,y  +  ^,z,f)- E,(x,y -^,z,t) 

Al  ■  Hm^(x,y,z)  L  2  2  2  2 

The  time!  derivative  operator  in  the  above  equation  is  typically  replaced  by  the 

finite!  difference  approximation  [Ref  5],  However,  just  like  in  1-D  and  2-D,  we  present  an 

alternate  approach,  such  that  the  approximation  of  the  field  time  variation  is  shown  to  be 

analogous  to  the  approximations  already  introduced  for  the  field  spatial  variation.  The 

above  equation  is  integrated  with  respect  to  time  to  get  an  approximate  equation  for  the 


magnetic  field  at  the  present  time  t 


(4.15) 


[£  Ey(x,y,  z-f,  x)dT  -  j'  Ey(x,y,  z  +  f,  x)dx  +  J'  Ez(x,y  +  f,z,  x)dx  -  j'  E,{x,y  -  f ,  z,  x)rfr] 

A  similar  integral  equation  can  be  written  for  the  magnetic  field  at  the  same  q)atial 


location  but  at  an  earlier  time  t-At 


Hx{x,y,  z,t-  At) »  - - ; - -X 

Al  •  Havgx(x,y,  z) 


(4.16) 


Ey(x,y,z-^,x)dx-^’J^  Ey{x,y,z+  ^,x)dx  +  \'J^  Ez{x,y-^-  ^,z,x)dx-\'^‘^  Ez(x,y-^,z,x)dx\^ 


Subtracting  the  two  equations  we  get 


ffxix,y,z,t)a  Hy,ix,y,t- At)  -  — - - - -x 

Al  ■  ilaygx(x,y,Z) 


(4.17) 


\^'^_^Ey{x,y,z-^,x)dx-\‘^_^  Ey(x,y,z+  f  ,'ry'r  + Ezix,y+  ^,z,x)dx-^'^_^Ez{x,y- ~,z,x)dt^ 


The  integrals  on  the  right-hand  side  can  not  be  evaluated  exactly,  because  the 
exact  tenq)oral  variation  of  the  electric  fields  within  the  At  interval  prior  to  t  is  grarerally 
not  known  (just  like  the  spatial  integral  for  the  magnetic  flux  could  not  have  been 
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evaluated  exactly  because  the  exact  variation  of  the  magnetic  field  within  the  contour  was 
not  known).  However,  the  integrals  can  be  evaluated  approximately  by  assuming  a  certain 
variation  of  the  electric  field  with  the  ten5)oral  variable  t.  The  sindplest  approach, 
consistent  with  the  assumptions  made  for  the  field  spatial  variation,  would  be  to  assume 
that  the  interval  At  is  small  enough  such  that  the  electric  field  can  be  assumed 
approximately  constant  within  At  and  ecjual  to  the  value  at  the  center  of  the  tune  interval 


(t-At,  t) 


Ey(x,y,z-f,x)'^ 

sEy(x,y,z-fj-f) 

for 

t-At<x<  t 

(4.18) 

Ey(x,y,z  +  j,x)f 

^Ey(x,y,z  +  jJ-f) 

for 

t-At<x<  t 

(4.19) 

Ezix,y  +  f,z,x)'^ 

^E^(x,y  +  f,z,t-j) 

for 

t-At<x  <  t 

(4.20) 

Ez(x,y-j,z,x)f 

-Ez(x,y-f,z,t-f) 

for 

t-At<x<t 

(4.21) 

The  above  represents  a  piece-wise  constant  approximation  of  the  electric  field 
variation  with  respect  to  the  temporal  variable  t.  The  approximate  expression  for  the 


magnetic  field  at  location  (x,y,z)  and  at  the  time  t  now  becomes 

HAx,y,  t)  »  >< 

The  above  equation  can  be  written  in  a  somewhat  different  form,  using  the  identity 

1  1  1  _  1 _ 1 _ 1 _ 1 _ L _ L  (AO'K 


1  1  1  _  1  11  1 _ L  =  v  J-J-  =  voTo— 


(4.23) 


where  the  free  space  velocity  of  propagation  is  denoted  V(„  the  intrinsic  impedance  of  free 
^ace  is  denoted  (Zo  «  3770),  and  daiotes  relative  permittivity.  Introducing  the 
grid  "propagation"  velocity 
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we  can  write  the  approrimate  "update"  equation  for  the  magnetic  field  as 
[Ey(x,y,z-f,t-f)-Ey(x,y,z  +  f,l-f)+E,(x,y+j,z,l-j)-Ey(x,y-j,z,l-^] 


(4.25) 


The  equation  simplifies  for  the  case  ofnon-magaetic  media  (ft  =1) 


H:zix,y,  z,  t)  «  Hx{x,y, z,  t-Ai)  -  Vo J  x 


(4.26) 


[Ey(x,y,  z-f.l-f)-E,(x.y,z+f,t-f)+E,(x,y+f,z,l-f)-E,(x,y-f,z,t-f)] 

The  update  equations  for  the  conq)onents  H^(x,y,z,t)  and  (x,y,z,t)  can  be 
derived  in  the  same  manner.  Since  the  only  differences  would  be  due  to  different  contour 
and  surface  orientations,  the  update  equations  for  the  components  Hy(x,y,z,t)  and 
(x,y,z,t)  can  be  obtained  from  the  update  equation  for  H,j(x,y,z)  by  cyclic  permitation  of 
indices.  The  contour  for  updating  Hy(x,y,z,t)  is  shown  below 


A  z 


Ez(x+dx/2,y,z,t) 


Ez(x-dl/2,yAt) 


Ex(x,y,z-dl/2,t) 


Figure  4.4.  A  Contour  for  Updating  B^. 


Following  the  same  procedure  shown  for  updating  con^onent,  or  singly  using 
the  cyclic  permutation  of  indices  we  get  the  update  equation  for  H^(x,y,z,t) 
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[e,(x  -  f,y,  z,  t  -  f )  -£r(x  +  ^,y,z,  t  -  f )  ■^E,(x,y,z  +  f ,  t  -  f)-E,(x,y,z-  f ,  t  -  f )] 
Finally,  the  contour  for  updating  H^(x,y,z,t)  is  shown  below 


Ex(x,y-dl/2^,t) 


Ey(x-dl/2,y,z,t) 


Hz(x,y,z,t) 


Ex(x,y-t-dl/2,z,t) 


Ey(x+dl/2,y,z,t) 


Figure  4.5.  A  Contoiu  for  Updating  H^. 

The  update  equation  for  H^(x,y,z,t)  is 

H^x,y,z,t)^H,(x,y,z.t-At)  - l'<.(lg;) >'  (««> 

[E.(x.y-f.z,t-f)-E.(x.y  +  f,z,t-f)+EAx*f.y,z.t-f)-E,(x-f,y,z,t-f)] 

In  most  cases  the  media  will  be  non-magnetic,  in  which  case  all  relative 
permittivity  averages  Avill  be  equal  to  1,  and  the  update  equations  will  simplify  accordingly. 

3.  Derivation  of  Electric  Field  Update  Equations 

We  will  start  with  the  Maxwell's  curl  equation  for  the  drculation  of  the  magnetic 
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H^(x,y,z,t)t  •  din  =  z(y,y,z^E:,(x,y,z,tyx  +Ey{x,y,z,t)y^  •  dsn  +(4.29) 

Ifsff  07 ]  • 

The  second  term  on  the  ri^t-hand-side  (the  conduction  current)  is  obviously 
going  to  cause  the  update  equations  for  the  electric  field  to  be  more  conq)hiicated  than  the 
magnetic  field  update  equations.  The  discrete  equivalent  of  the  magnetic  field  circulation 
needs  to  be  determined  for  three  sample  contours  a  contour  in  a  plane  parallel  to  the 
yz-coordinate  plane,  a  contour  in  a  plane  parallel  to  xz-coordinate  plane,  and  a  contour  in 
a  plane  parallel  to  the  xy-plane.  The  first  contom  provides  an  update  equation  for  the 
x-directed  conq)onent  of  the  electric  field,  the  second  contour  provides  an  update  equation 
for  the  y-directed  component  of  the  electric  field,  and  the  third  contour  provides  an  update 
for  the  z-directed  component  of  the  electric  field.  Again,  it  is  sufficient  to  derive  the 
update  equation  for  one  electric  field  component  only  (say  the  x-conponent)  since  the 
other  equations  can  be  obtained  by  a  simple  cycHc  permutation  of  indices. 

"  A 

i 

Hz(x,y-(fy/2,z,t) 

X 

^  Hy(x,y,z-dz/2,t) 

Figure  4.6.  A  Contour  for  Updating  E^. 
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A  local  coordiaate  system  (\(/,  <;)  will  be  used,  with  the  origin  at  the  center  of  the  contour. 
Any  point  (y’,z')  within  or  on  the  contour  Cjj(or,  equivalently,  any  point  on  the  surface  Sjj) 

can  be  specified  by  its  local  coordinates 

y  =y  +  \]f  z^=z  +  q  where  -^<q<+^,md  -^<v|/<+^ 

The  circulation  of  the  magnetic  field  around  ,  assuming  counter-clockwise 

reference  direction  such  that  the  normal  to  the  surface  Sjj  is  in  the  direction  of  the  electric 
field  Ex(x,y,z,  t),  is  approximately 

^Hy(x',y,z',  t)y  +Hz{x'y,z',  ffz  ]  •  «  (4.30) 

[Hy(x,y,z-j,t)-Hy(x,y,z  +  f,t)+H^(x,y  +  f,z,t)-Hz(x,y-f,z,t)]-Al 
The  electric  flux  through  the  surface  Sjj  can  only  be  evaluated  approximately,  snice  the 

exact  way  that  the  electric  flux  density 

Dxix,y,  z,  t)  =  s(x,y,  z)Ex(x,y,  z,  t)  (4.3 1) 

varies  with  y  and  z  is  not  known  a  priori.  First  we  re-write  the  flux  integrals 

e(x^,/,  z',  t)Ex(x',y, z',  t)  x  *dsH=  (4.32) 

\A  J  J  +  V. -Z  +  <i)Ex{x,y  +  vp, z + <;,  0  X  •  X  dxydq 

2  2 

Ssff  •dsH=  (4.33) 

^4^  r+~  — ^  ^ 

JJJJ  o(x,y  +  \i/,z  +  q)Ex(x,y  +  \\/,z+q,t)x  •  xd^idq 

2  2 

Next,  we  assume  that  the  contour  side  A1  is  small  enough  such  that  the  electric 
field  E(x,y+v}/,z+^,t)  within  the  contour  may  be  assumed  constant  and  equal  to  the  value  at 
the  contour's  center  E(x,y,z,t).  This  assunqrtion  yields  approximate  expressions  for  the 
flux  integrals 

s(x,y  +  \y,z  +  q)Ex(x,y  +  vp, z  +  q, t)t  (4. 34) 

— >  r+— 

•dsff^  Ex(x,y,  z,t)-j_lj  J  e(x,y  +  vp, z  +  q)d\ydq 

2  2 
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o(x,y  +  \\i,z+(^Ex{x,y  +  \^,z+c„t)x  (4.35) 

•dsH^ Ey{x,y, _l  <y(x,y + h/, z + q)d\ifdq 

2  2 

or,  using  the  permittivdty  and  the  conductivity  averaged  in  the  y  and  z-directions 

8(x,y  +  vj/, z  +  <^Ex(x,y  +  vj/, z  +  ?, t)x  •dsn^  {Mf  •  £»(x,7,z,  t)  •  Za.gx{x,y,z)  (4.36) 
a(x,y  +  M/,  z  +  <;)i:;.(x,y  +  ^,zyq,  t)x  •dsn^  {Mf  •  Ex{x,y,z,  t)  •aavgx(x,y,z)  (4.37) 

The  second  curl  equation  of  Maxwell  can  now  be  replaced  by  an  approximate  equation 

[/7/x,y,z-  f)-i^^(x,y,z  +  f)  +Hz(x,y  +  j,z,t)-H^(x,y-j,z,t)]  -  A/ «(4.38) 

^{[Savgx(x,y,z)£:^(x,y,z,0](A0^}  +(A/)^  •<t.n«;r(x,j,z)  •Ex(x,y,z,t) 

The  approximate  curl  equation  can  be  re-written  as 

§;{E.(x,y,z,0)  (4.39) 

- L - -[//,(x,7,z  -  f ,  0 -Hy(x,y,z+j,  0  +Hr(x,y+j,z,  i) -H4x,y-f,z,  ()] 

A/  *  Zavgx\P^,yy^) 

i  ^avgx(y,y,Z^ 

-.^(x,y,z) 

The  above  equation  is  integrated  with  respect  to  time  to  get  an  approximate  equation  for 
the  electric  field  at  the  present  time  t 


E:z(x,y,z,t) 


(4.40) 


A/  •  ^avgjci^yyy 

[I'  ffy(x,y,  2  -  f,  't)rfT  -  £  ffy(x,y,  z  +  f,  x)dx  +  J'  Hz{x,y  +  f  ,z,  x)dT  -  £  Hz{x,y  -  f ,  z,  x)rfx] 

_a^.(x,y,z)  » 

Za^{x,y,z)  Jo 

Similar  equation  can  be  written  for  the  electric  field  at  the  same  spatial  location  but 


at  an  earlier  time  t-At 


Ex(x,y,  z,t-  At)  «  — - - -  X 

AI-Eaypcix,y,z) 


(4.41) 


HyiXyyyZ-^yXyix-^’^^^  Hy{x,y,z  +  ^,x)dz  +  \‘J^  Hz(x,y+^,z,x)-\‘'^  H:(x,y- ^,z,x)dx^ 

Oav^{x,y,z)  f/-A/ 
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Subtracting  the  two  equations  we  get 


(4.42) 


\J^H,(>c,y,z-^,zyh-\^_^H^x.y.z+^.z)dzy^^H.Xx.y*^,z,z)<h-\H.(x,y-^,z,%)ii'\ 

-  ’z^(x,y,z)  r,  ^  ^  y 

&a^{x,y,z)  h-u 

Assuming  that  the  interval  At  is  small  enough  such  that  the  magnetic  fields  can  be 
assumed  approximately  constant  within  At  and  equal  to  the  value  at  the  center  of  the  time 


interval  At 


/f,(W  +  f,T)»ft(w  +  f>'-f)  for  '-A<StS( 

H,(x,y,z-f,z)«H,(x,y,z-f,t-f)  for  l-Al<z<t 
H,(x,y+f,z,x)r>H^x,y-f,z,l-f)  for  l-AI<x<t 
H,(x,y-f,z,z)<^H,(x,y-f,z,l-f)  for  t-Atizit 


we  get 


Ex(x,y,  z,  t)  «  Ejcix,y,  z,  t-At)  + 


A/  •  ZavgxiX,y,z) 


(4.43) 

(4.44) 

(4.45) 

(4.46) 


(4.47) 


[Hy(x,y,  z  -  f ,  t  -  f )  -  Hy{x,y,z  +  +Hx(x,y  +  f  ,z,  t-  f )  -  Hz(x,y  f )] 

_Gaygx(X,y,z)  « 

^avgx(x,y,z)  t  &t 

The  last  integral  term  can  be  expressed  in  different  form  as  shown  below 


The  above  equation  can  be  written  in  more  compact  manner  by  using  the  identity 

1  1  1  1  1 _ (4  49) 


i  —  J_-L  -  ^ _ X  1  _  X  ^  1  _  y  1 

^~^oer- yp^VSoS.  0  0s. 


vsdiere  the  firee  space  velocity  of  propagation  is  denoted  v„,  the  intrinsic  impedance  of  fi-ee 
space  is  denoted  Zg  (Zo  «  377Q)  and  s^  denotes  relative  permittivity.  Introducing  the  grid 
"propagation"  velocity 


V 

Vgrid  - 


(4.50) 


104 


we  can  write  the  approximate  equation  for  the  electric  field  "updates"  as 

E,{x,y,z.,)^E.(x.y,z,x)  (4.51) 


The  average  of  the  electric  field  between  t-At  and  t,  denoted  can  not  be 

calculated  exactly  because  the  exact  tenq)oral  variation  of  the  electric  field  within  the  At 
interval  is  not  known  a  priori  However,  assuming  a  linear  variation  within  At  interval,  the 
terq)oral  average  of  the  field  is  the  average  of  the  field  values  at  t-At  and  t 

(4.52) 

Substituting  the  time-average  Ey  ^(^,)(x,y,t)  into  the  update  equation  we  get 

Eyc(x,y,z,  t)  «  E:,(x,y,z,  t  -  At)  +2o^; - T  x  (4.53) 

''gnd  Sr^x{X,y,Z) 

[Hy(x,y,  z  -  J,  ^  -  f )  -  /r^(x,y,  X  +  f,  t  -  f)  +Hz{x,y  +  f,z,t-  At)  -  H^(x,y  -  f,  z,  ^  -  f )] 

^  Vo  1  _  E:cix,y,z,t-At)+Ex(x,y,z,t) 

-^»V^e,^(V,i)  ■  - 2 - 

which  give  the  update  equation  for  the  x-cocDqponent  of  the  electric  field  as 

^  1  Vo  ZQ  Al-Oayg::{x,y,z) 

Ex(x,y,  z,  t)  « - ^  t  -  A/)  +  (4.54) 

J  ^  1  Vq  Zq  •  a/  •  Ogvgxix^y,  z) 

2  ^r/ivgx(x,y,z) 

Z  ^0  1 

_ V^”<(  ^r,aygx(x,y,z) _ ^ 

1  ^  1  Vq  Zo  Al-aa^(x,y,z) 

'^2'^ grid  e,.^(x,y,r) 

[Hy(x,y,z-f,t-f)-Hyix,y,x  +  f,t-^)+H,(x,y  +  f,z,t-At)-Hz(x,y-f,z,t-f)] 

The  equation  simplifies  greatly  for  the  case  ofnon-conductive  media  (o=0) 
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Ex(x,y,z,  t)  « Ex(x,y,z,  t-  At)  +Zo 


(4.55) 


^^^vgnds,^^,(x,y,z)  ^ 

\Hy{x,y,  z-f,t-f)-  Hy(x,y,  z  +  f,t-f)  +H,(x,y  +  f,  z,  /  -  At)  -  H,(x,y  -f,z,t-f)] 
In  the  case  of  free-space  (8=1),  the  equation  simplifies  fiirther 


Vn 

Ex(x,y,z,  t)  «  E:c(x,y,z,  t  -  At)  +2q—  x 


(4.56) 


[Hy(x,y,  z  -  f ,  t  -  f )  -  Hy(x,y, x  +  f ,  /  -  f )  +  H,(x,y  +  f,z,t-  At)  -  H,(x,y  -f,z,t-  f )] 
The  update  equations  for  the  other  two  electric  field  conq)onents  Ey(x,y,z:,t)  and 

(x,y,z,t)  can  be  derived  in  the  same  manner.  Since  the  only  differenc'es  are  in  the  contour 
and  surface  orientation,  the  update  equations  for  Ey(x,y,z,t)  and  E^  (x,y,z,t)  can  be 
determined  by  a  simple  cyclic  permutation  of  indices,  as  shown  on  the  next  page.  The 
contour  for  updating  Ey(x,y,z,t)  is  shown  below. 


Hz(x+dx/2,y,z,t) 


Hz(x-dl/2,yAt) 


Hx(x,y^-dl/2,t) 


Figure  4.7.  A  Contour  for  Updating  Ey. 

Following  the  same  procedure  as  shown  for  the  E^(x,y,z,t)  update  equation,  or  by 
a  sirople  cyclic  permutation  of  indices  we  obtain 
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Ey(x,y,z,t) 


J  1  Vq  Zq  •  dkl  •  Oavgy(x,y,Z) 

+  (4.57) 

I'^grid  er.a^gy(x,y,z) 


7  ^  _ j- 

*^Vg77rfSr,avfy(y,J^,^)  ^ 

j  ^  1  Vq  Zq  •  •  Oavgy(x,y,Z) 

Ingrid  er,avgy(X,y,z) 

[H,ix  -  f  ,7, 2,  t-M)-  Hz(x + f,y,z,  t-j)+  Hx(x,y,z  +  -H:,(x,y,  z-^,t-f)\ 

The  equation  simplifies  greatly  for  the  case  of  non-conductive  media  (a=0) 


Vo  1 

EAx,y.z.t)  »  Ey(x,y,z,t-Ai) 


(4.58) 


[Hz(x  -  J,y,  z,  t-At)-  H,(x  +  J,y,  z,  ?  -  f )  +  Hxix,y,z + f,  /  -  f )  -  H:c(x,y,  ^  -  f ,  ^  -  f)] 
In  the  case  of  firee-space  (s=l),  the  equation  sinq)Iifies  fiirther 


Ey(x,y,z,l)s! Ey(x,y,z,t-Al)  yZo-^x 


(4.59) 


— y  Vgrid  ^ ' 

[h,(x  -  f,y,  z,  t-At)-  Hz{x  +  f  z,  t  -  f )  +  H,c(x,y,z + f,  t  -  f )  -  Hx(x,y,  z  -  f,  t  -  f )] 
The  contour  for  updating  E^(x,y,z,t)  is  shown  below. 


Figure  4.8.  A  Contour  for  Updating 


Following  the  same  procedure  as  shown  for  the  Ex(x,y,z,t)  update  equation,  or  by 
a  siinple  cychc  permutation  of  indices  we  obtain 

j  1  Vo  ^0  '  '  ^avgz(x,y,Z) 

(4.60) 

'^2'^ grid  Sr^vgz(x,y,z) 

2  ^0  1 

”v^r,JSr.^Cy,y,z)  _ 

j  ^  1  Vo  ^0  '  ^ avg:(x,y,Z^ 

2^  grid  Sr^avgz(x,y,Z^ 

[H,{x,y  -  f ,  z,  ?  -  f )  -  H,{x,y  +  f ,  z,  f  -  f )  +Hy{x  +  f  ,y,  Hy(x  -  f,y,z,  t  -  f )] 

The  equation  simplifies  greatly  for  the  case  of  non-conductive  media  (cr=0) 

E,(x,y,z,i)«E,(x,y,z,t-M)  (4.61) 

[H.(x,y  -  f  ,z,  /  -  f )  -  ft(x,y + f  ,z.  /  -  f )  +(?,(x  +  ^,y,  z,t-f)-  Hy(x  -  f,y,z,  ( -  f )] 
In  the  case  of  firee-space  (s  =1),  the  equation  anq)lifies  further 


Ez(x,y,z,  t)  «  Ez(x,y,z,t-M)  +Zo^x 


4.62) 


[Hy(x,y-^,z,t-f)-Hy(x,y  +  f,z,t-f)+Hy(x+f,y,z,t-f)-H,(x-f,y,z,t-f)] 
B.  TRANSPARENT  GRID  TERMINATION  IN  3-D 


1.  3-D  Grid 

The  electric  and  magnetic  field  update  equations  have  been  derived  using  local 
coordinate  systems,  with  origins  at  the  centers  of  the  magnetic  and  electric  contours, 
respectively.  These  equations  now  need  to  be  "converted"  to  a  global  coordinate  ^stem, 
that  is  to  a  3-D  grid  of  equi-distant  sampling  points.  We  assume  that  our  domain  is  an  L 
by  L  by  L  cube  discretized  using  a  uniform  grid,  that  is  that  the  discretization  step  A1=L/N 
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is  not  nhanging  with  position.  The  electric  and  magnetic  field  san:q)ling  points  (nodes)  are 
"interleaved",  that  is  offeet  by  Ah2  from  each  other  in  x,  y,  and  z-directions.  The  electric 
field  con:q)onents  are  located  along  the  sides  of  "electric"  cells  while  the  magnetic  field 
components  are  located  along  the  sides  of  "magnetic"  cells. 


The  electric  cells  are  the  cells  of  constant  e  and  o,  while  the  magnetic  cells  are  the 
cells  of  constant  |j..  The  permitivity  and  conductivity  can  vary  from  one  electric  cell  to 
another  in  arbitrary  maimer.  The  permeabilty  can  vary  from  one  magnetic  cell  to  another  in 
an  arbitrary  manner  as  well.  The  electric  cells  are  shifted  by  hM2  in  all  three  coordinate 
directions  with  respect  to  corresponding  magnetic  cells.  The  placement  of  field 
components  along  the  cell's  sides  automatically  provides  for  the  continuity  of  tangential 
components  of  electric  fiels  across  media  interfaces,  while  the  spatial  shift  of  electric  vs. 
magnetic  cells  provides  for  electric  and  magnetic  field  linkage  (coupling).  Domain 
boundaries  can  coincide  with  either  the  fiices  of  electric  cells  or  the  faces  of  magnetic  cells, 
but  not  both  (because  of  the  Al/2  spatial  shift  between  the  two).  Although,  in  principle, 
either  electric  or  magnetic  cells  can  be  alligned  with  the  grid  boundaris,  we  will  select 
electric  cell  alligment  with  the  domain  boundaries. 
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Figure  4.10.  Yee  Electric  Cell  and  Associated  Field  Con5)onents. 

A  grid  of  (N  by  N  by  N)  electric  cells  has  (NXN)(N)  E^,  E^,  and  E^  electric  field 
nodes,  (N-lXN-lXN-2)  Hy,  and  magaetic  field  nodes  (Just  the  opporite  would  hold 
for  alhgning  magnetic  cells  with  domain  boundaries).  3-D  electric  and  magnetic  field 
"update"  equations  can  now  be  "converted"  to  electric  and  magnetic  field  grid  "iq>date" 
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equations  by  replacing  the  variables  x,y,z  and  t  with  the  grid  coordinates  of  the  electric 
field  spatial  and  temporal  sairq)ling  points 

X  ->•  /A/,  i  =  0, 1...N  y  -» jM,  j  =  0, 1...N  z  M/,  A:  =  0,  l..A^  t  ->  wA/,  k  =  0, 

The  notation  can  be  sinaplified  by  omitting  the  common  A1  and  At  terms,  and  u^g 

a  superscript  for  the  index  of  the  temporal  sanq)ling  point.  The  E  field  grid  equations  for 

conductive  media  are  obtained  from  equations  (4.54, 4.57, 4.60) 

1  'Z^.A/.rr  (i  i 


■■ 


j  _  1  Vq  Zq  •  Al  •  Oavgxijfjt 

^  ^  1  Vo  Zo-Al-G^^(iJ,kr^  ^ 

grid  £  rfivgx  (j,j, 


^  Vo  1 

^0,.  . - r 


^ grid  Srfivgx(j,j\  ^ 

^  ^  1  Vq  ZQ-M-<5avgxiiJ,k) 

2,^  grid  £^^,^(x,y,z) 

H7hiJ,k-\)-H7\iJ,k+\)  +H7kiJ +  hk)-H7^\iJ -\,k 

j  _  1  Vq  ■^O  '  ‘  ^avgy(i,J,k) 

J  I  1  ^0  ^0  ‘  '^oygyv^J^fQ 

2  E r,avgy(i ,7,  A) 

^  Vo  1 
Zot: — : - 7r":  “ 


(4.63) 


(4.64) 


_ S  r,avgy(jyjy  ^ 

I  ^  1  Vq  Zq  •  Al  •  Oavgyjiyjy 

2  ^grid  g r^gyiij^ 
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(4.65) 


\  _  1  VQ  Zq  •  AI  •  OgygzQj^ 

I  I  1  ^0  -^0  ’  ‘  ^avgz\}J^  k) 

Ingrid  Zr,avgz(i,J,K) 


^g’dd  Zr,avgz(j^jt^  ^ 

J  _l_  1  Vq  Zq  •  Al  •  Gavgz(i,j,k^ 

2^gnd  Sr,m>gziiJ,k) 

H7^(iJ-\,k)-Hx  ^{iJ  +  \,k)+H7\i  +  \j,k)-Hy  \i-\j,k) 
The  equations  surq)lify  greatly  for  the  case  of  non-conductive  media  (o=0) 


E"^(i,j,k)aE"^\i,j,k)  +  Zo 


^grid  £r,avgxij,j,  k) 


H7khj,k-\)-Hf\iJ,k+  \)  +Hf\iJ  +  \,k)-H7hiJ -  hk) 


Eyi^J,k)»E"^{iJ,k)  +  Zo 


^grid  £r^gy(i,J,  k) 


H”z\i-  \j,  k)-Hz  ^  (/•  +  U,  k)  +ifr^ (ij,  k+\)-  h7^  QJ,  k-\) 


(4.66) 


(4.67) 


E”{iJ,k)»ET\iJ,k) 

Hx  \iJ-\,k)-Hx  \iJ  +  \,k)+H7''{i  +  \j,k)-H7^ii-\j,k) 

In  the  case  of  free-space  (s  =1),  the  equations  sircplify  further 

E%i,J,k)^Er'(t,j,K)  +Zo;^x 


E",(iJ,k)«E^'(iJ,k)  +Zov^x 

j,  k)-H:'’  (/  +  U  k)+HT'  (t,j,  * + i)  -  (a  k-\) 


E"M,k)^ErKi,J,k) 


Hf\i,j~\,k)-H7'(i,i  +  \,k')+Hf\i*\j,k)-H7'\‘-\M 


(4.68) 


(4.69) 


(4.70) 


(4.71) 
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Similarly,  the  H-field  update  grid  equations  are 


ETku,  k-\)-  eT'  (ij,  k  +  ^)  +Ef\iJ  +  h  k) -ET\iJ  -  \,k) 


m(iJ,k)^H;-\iJ,k)  -  Yo 


Vo  1 

n,(/,y.  4)' 


fip  (/ -\,J.k)-ET\l  +  \,J,k)*  eP  (U  4 + i)  -  eT%J,  *  -  5) 


(4.72) 


(4.73) 


Ep(i,j-\,k)-Ep(i,J*\>'^)+E",  \i  +  \j,k)-E",  ^(t-\j,k) 


(4.74) 


2.  3-D  Grid  Termination 


The  grid  equations  derived  previously  are  valid  for  all  the  nodes  except  for  the 
nodes  on  the  domain  boundary.  The  reason  is  that  the  field  nodes  on  the  boundary  have  at 
most  five  nearest  neighbor  nodes  while  the  nodes  within  the  domain  boundaries  have  six 
nearest-neighbor  nodes. 
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because  the  field  couponent  on  the  boundary  shown  in  the  figure  above  can  not  be 
updated  using  the  grid  equations  discussed  so  &r.  Extending  the  ideas  of  transparent  grid 
termination  (TGT)  the  discrete  boundary  impulse  response  (DBER.)  applied  for  1-D  and 
2-D  problems.  Similarily  to  the  2-D  TGT  concept,  the  update  equations  on  the  botmdary 
can  be  devised  based  upon  convolutions  of  the  fields  on  the  just  inside  boundary  and  the 
pre-calculated  impulse  responses  from  the  just  inside  boundary  to  the  boundary.  The 
fields  on  the  boundary  will  be  e?q)ressed  as  a  superposition  of  rehouses,  as  illustrated  in 
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the  figure  below  for  top  plane  of  the  boundary.  The  same  applies  to  other  plane  of 
boundary  as  well. 


the  node  on  the  boundary  has 


Figure  4.12.  Fields  on  the  Boxmdaiy  as  a  Superposition  of  Re^onses. 

The  superposition  can  be  also  modelled  using  the  concept  of  a  multiport.  The 
nodes  on  the  grid  boundary  have  three  field  components  and  those  field  components  can 
be  conadered  as  output  ports  of  a  linear  multi-port.  Assuming  that  the  grid  is  a  cubic  with 
(NXNXN)  nodes,  the  total  number  of  output  ports  will  be  3[2(NXN)+(N-2)4(N-1)].  The 
field  components  on  the  nodes  just  inside  the  grid  boundary  will  be  considered  as  input 
ports  of  the  linear  multiport.  There  will  be  a  total  of  3[2(N-2)(N-2)+(N-4)4(N-3)]  input 
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ports.  The  equivalent  multiport  model  in  2-D  TGT  can  be  used  in  the  same  manner.  The 
input  and  output  fields  can  be  either  electric  or  magntic  fields.  Furthermore,  the  input  and 
output  fields  need  not  be  of  the  same  "kind",  that  is  the  input  fields  could  be,  in  principle. 
E-fields  and  the  output  fields  could  be  H-fields,  or  vice  versa.  How^ever,  we  will  select  the 
input  and  output  fields  to  be  of  the  same  "kind",  that  is  either  all  E  or  all  H  fields,  such 
that  the  impulse  responses  h^(t)  will  be  dimensionless.  The  equations  will  be  derived  for 
E  field  model,  that  is  with  electric  fields  on  the  grid  boundary.  Dual  equations  for  the 
magnetic  fields  model  can  be  obtained  shnply  by  replacing  E's  with  ITs.  There  are  six 
planes  for  the  cubic  grid  boundary.  However,  by 
maldng  an  arrangement  for  the  input  and  output  in 
TGT  model  as  shown  below. 


splitting  the  cubic  into  N  pages  and 
the  N  pages,  we  can  simplify  the  3-D 
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Figure  4.13.  Page  Arrangement  for  3-D  Cubic  Cell 

The  boundary  field  coirqponents  for  the  first  page  have  the  electric  field  E^(iJ,l,tX 
Ey(ij,l,t)  and  E/ij,l,t).  The  boundary  field  conqjonents  for  the  second  page  through 
(N-l)th  page  have  the  electric  field  components  like  this.  First,  for  the  x  components,  the 
2nd  Page  has  [E,^(lj,2,t),  E,(i,N,2,t),  E,(Nj,2,t),  E,(i,l,2,t)],  nth  (2<n<N-l)  page  has 
[E,.(lj,n,t),  E,(i,N,n,t),  EiNj,n,t),  E,,(i,l,n,t)]  and  (N-l)th  page  has  |E,(lj,N-l,t), 
E,,(i,N,N-l,t),  E,,(Nj,N-l,t),  E,,(i,l,N-l,t)].  For  anq)Hcity,  the  y,z  components  have  the 
gamft  global  coordinates  used  in  the  x  corcponoats.  Finally,  the  boundary  field  con:q)onents 
on  the  Nth  page  have  E,,(ij,N,t),  Ey(iJ,N,t)  and  E/ij,N,t).  The  limits  on  the  values  of  ij,k 
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on  the  boundary  are  1  and  N.  Now,  the  field  components  on  the  nodes  just  inside  the 
boundary  exist  in  the  second  page  throu^  the  (N-l)th  page.  The  second  page  has 
Ey(ij,2,t)  and  E^(ij,2,t).  For  the  x  coi^poonents,  the  third  page  has  [E,j(2j,3,t), 
E,(i,N-l,3,t),  E,(N-lj,3,t),  E,(i,2,3,t)],  mth  (3<m<N-2)  page  has  [E,(2j,m,t), 
E,^(i,N-l,m,t),  E,^(N-lj,m,t),  E,^(i,2,m,t)].  The  y,z  components  have  the  same  global 
coordinates.  Finally,the  (N-l)th  page  has  E^(ij,N-l,t),  Ey(ij,N-l,t)  and  E/ij,N-l,t). 
However,  the  limits  on  the  values  of  ij,k  on  the  nodes  just  inside  the  boimdary  are  2  and 
N-1.  Because  the  each  node  on  the  boundary  has  three  field  components,  the  electric  field 


vectors  for  the  grid  on  the  boundary  can  be  formed  as  below 

Ex^boundaryif)  “  \Eix_\st jfage  Exjlnd _page  ••• 

•  •  ExJJ^-Vyh Jfage  ExJ^th j>age\^ 

(4.75) 

Ey_Jboundary(J)  ~  \Ey_\st _page  Eyjind jfage  •- 

Eyjj^-iyh Jfage  Eyjslth jfage\^ 

(4.76) 

Ez^boundaryif)  ~  [Ez^lst jfage  Ezjlnd jfage  •• 

*•*  EzJ^N-iyth Jfage  EzJsJrh jfage] 

(4.77) 

where  each  page  field  element  forms  a  row  vector  with  the  fields  at  the  nodes  on  the 
boundary  and  T  indicates  transpoation  (the  vector  is  a  colxmm  vector,  but  it  has  been 
written  as  a  transpose  of  a  row  vector  to  save  ^ace).  In  an  analogous  manner  an  electric 
field  vector  can  be  formed  for  the  nodes  just  inside  the  grid  boundary 


Ex justjmideif)  ~  [Exjlnd jfage •••• 

....ExJN-V^h Jfage] 

(4.78) 

EyJustjnsideif)  ~  [EyjZnd jfage •••• 

...£y_{fj^\)th Jfage] 

(479) 

Ez Justjnsideif)  ~  \EzJlnd jfage . 

(4.80) 

Ejustjnsideif)  “  \Ex justjmideif)  Eyjust_ 

imideif)  Ez JustJnsideif)] 

(481) 

vdiere  each  page  field  element  forms  a  row  vector  with  the  fields  at  the  nodes  just  inade 
the  boundary.  The  electric  field  on  the  grid  boimdary  Eb„j,^(t)  can  be  expressed  as  a 
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convolution  of  the  electric  fields  at  the  layer  just  inside  the  grid  boundary  i^ij^(t)  and  a 

matrix  of  ii^pulse  re^onses  H(t). 

ExJjoundaryif)  ~  f^x(t)  *  Ejustjnsideif)  (4, 82) 

Eyjyoundaryif)  Elyif)  *  Ejustjnsideif)  (4-83) 

E z_bounfiaryif)  *  Ejustjnsideif)  (4-  84) 

wheie.  the  2(N)(N)^(N-2>im-l)  by  3[2(N-2)(N-2)+(N-4)4(N-3)]  matrix  of  impulse 


Hy(t)  and  H^(t)  have  the  same  matrix  size  as  E^(t).  The  electric  field(Ex,  Ey,  Ez)  at 
a  boundary  node  (identified  by  a  particular  value  of  the  index  m)  can  be  written  as  the 
product  of  the  m-th  row  of  the  impulse  response  matrix  H(t)  (Hx(t),  Eij,(t),  H^(t))  and  the 
column  vector  of  the  electric  fields  just  inside  the  boundary  (indentified  by  index  n) 

Ex_boundary(jn,  t)  =  hx  * Ejustjnside(a,  f)  =  (4.86) 

2  n=l  Jq  Ejustjnsideif}  t)  *  ~ 

EyJ,oundatyOf,t)=^„-l  hy  *  Ejusrjnside(n,i)  —  (4-87) 

S  „=i  Jq  Ejustjnside  (f,  t')  ’  msi(t  ~ 
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Ez_boundaryij^,^  —  ^Jn=\  *EjustJnside(n,f)  =  (4.88) 

2  Jq  Ejust_inside(f^}  t)  "  mftit  ~  T)^^T 

The  discretized  forms  of  the  above  equations,  involving  samples  of  electric  fields  at  t=kAt, 
are  the  sums  of  discrete  convolutions.  Since  a  discrete  convolution  is  itself  a  sum,  the 
result  for  the  boundary  nodes  will  involve  double  summation. 

=  s:3'  2^  ■  t&t,  (4-89) 

Ey^bouHdaiyi'^)  ~  2^  Efj„!t_haidA^  '  hymjl  (4.90) 

=  2:5'  2^  -hiz,-  (4.91) 

The  convolution  equations  express  the  fields  at  the  edges  as  weighted  sums  (the 

wei^ting  coefi&caits  are  the  values  of  the  discrete  boundary  hrpulse  req)onses  h^)  of  the 

time  histories  of  the  fields  "just  inside"  the  grid.  In  that  respect  the  concept  and 

procedure  used  in  3-D  TGT  is  the  same  as  those  used  in  2-D  TGT.  Also,  the  impulse 

responses  h^(t)  converge  to  zero  relatively  fast  v^Mch  effectively  reduces  the  number  of 

significant  terms  in  the  convolution  sums.  The  convolution  sums  can  thus  be  truncated 

such  that  only  a  certain  number  of  the  most  recent  values  of  the  electric  fields  just  inade 

the  boundary  are  needed.  The  discrete  boundary  inq)ulse  req)onses  hj^(kAt)  need  to  be 

determined  (and  saved)  only  once,  just  like  in  1-D,  2-D,  for  a  selected  grid  velocity 

Vg^j=Al/At.  The  impulse  responses  can  be  stored  in  a  rectangular  matrix  form.  A 

particularly  suitable  matrix  form  would  have  [2(N)(N)+^-2)4(N-l)]  rows  and 

3[2(N-2)(N-2)+(N-4)4(N-3)]Nb  columns(  because  the  field  components  at  the  nodes  just 

inside  the  boundary  consist  of  Ex,  Ey  Ez  and  the  nmnber  of  each  elctric  field  is 
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2(N-2XN-2)+(N-4)4(N-3)])  ,  where  denotes  the  "truncated"  length  of  the  discrete 
boundary  mq>ulse  responses.  This  form  is  suitable  because  the  boundary  electric  fields  can 
then  be  calculated  as  a  product  of  such  a  matrix  and  a  column  vector  constructed  of 
3[2(N-2XN-2)+(N-4)4(N-3)]  subvectors  of  length  N^.  The  subvectors  represent  the 
most  recent  values  of  the  electric  field  at  nodes  just  inside  the  boundary.  The  choice  of 
depends  on  the  desired  amplitude  "resolution".  These  subvectors  are  updated  at  each  time 
step  by  shifting  all  of  their  values  down,  that  is,  "fiither  into  the  past",  by  one,  discarding 
the  last  value,  and  updating  the  first  subvector  element  with  the  most  recent  field  value 
calculated  via  standard  grid  update  equations.  The  double  summations  are  thus  efficient^ 
executed  as  matrix-vector  multiphcations. 

In  the  3-D  TGT,  the  important  thing  is  how  to  determine  the  diso-etized  impulse 
responses.  The  discretized  impulse  responses  are  determined  in  a  manner  analogous 
to  the  1-D,  2-D  cases,  using  the  discrete  equivalent  of  the  Dirac  deha  function  \\diich  we 
will  denote  as  5^  In  order  to  determine  the  discrete  boundary  impulse  response  for  a 
particular  input  port  n,  all  input  ports  except  the  n-th  input  port  are  set  to  zero  and  the 
Dirac  impulse  is  applied  to  the  m-th  input  port.  Recording  all  the  outputs  from  t=0  to 

t=N^At  provides  us  with  2(NXN)+(N-2)4(N-1)  discrete  boundary  impulses  of  duration  N^. 

Eljustjr.iae(n)  =  h^  (4.92) 

Eljustjnsideip)^^  for  p  =  -1^ +iN -  -3)  and  56 /I  (4.93) 

E';justjnsiae(p)  =  0,Elj„^j„,,,,(p)  =  0  for  /;=  1...2(Ar-2)2 +(iy-4)4(Ar-3)  (4.94) 

^0  determined  in  the  same  manner  as  the 

^^j„,OWe(«)  =  5*  (4.95) 
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Eyju^j,^ide(p)  =  ^  for  p  =  \...2iN-iy+(N-A)A{N-3)  and  p^n  {AM) 
E].j.sUnsiM  =  ^Xj.ur>^iM  =  ^  for  p  =  l..2{N-2y  +{N-mN-2){A.91) 

El ju^jnsidein)  =  6*  (4.98) 

Eljustjnside(p)  =  ^  for  p  =  l...l{N -if +{N - A)A{N -3)  and  p^n  (4.99) 
EljustJnside(P)  =  ^,Elj„sUr^,M^  fOT  p  =  -if  ^{N  -  A)A{N -3) 

(4.100) 

The  each  process  is  repeated  2(N-2XN-2)+(N-4)4(N-3)  times  and  the  results  for 
the  each  case  are  stored  in  the  2(N)(N)+(N-2)4(N-1)  by  3[2(N-2)(N-2)+(N-4)4(N-3)K  ( 
[E^  Ey  EJ  field  vector  just  inside  the  boundary)  matrix.  This  matrix  is  then  used  to 
implement  the  transparent  grid  termination  via  the  matrix  multiplication  by  the 
3[2(N-2)(N-2)+(N-4)4(N-3)]N^  column  vector  of  the  time-histories  of  the  fields  just  inside 
the  boundary.  Depending  on  the  values  of  and  N,  there  will  be  input  nodes  that  are 
sufficiently  far  from  an  ouq)ut  node  such  that  the  time  it  takes  for  an  input  in^ulse  to 
propagate  to  the  ouput  port  exceeds  N^.  In  such  a  case,  the  contribution  of  such  an  input 
node  to  the  output  node  is  known  to  be  zero  in  advance.  This  may  be  used  to  reduce  the 
number  of  operations  in  TGT  implementation.  The  process  of  obtaining  the  3-D  discrete 
boimdary  inq)ulse  responses  is  carried  out  on  a  "large"  grid  whose  size  should  be  at  least 
(N+N^  by  N+N^  by  N+NJ,  or  equivalently,  the  distance  from  the  TGT  boundary  to  the 
boundary  of  the"large"  grid  should  be  at  least  NbAy2.  The  termination  of  this  "larger"  grid 
(in  typical  applications  N»N^  and  the  "larger"  grid  would  be  only  incrementaly  larger)  is 
immaterial,  since  the  reflections  off  its  boundary  do  not  arrive  at  the  output  nodes  (where 
they  would  have  corrupted  the  impulse  responses)  before  the  time- stepping  has  been 
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terminated.  When  calculating  the  discrete  boxmdaty  impulse  responses  (DBIR),  the  nodes 
interior  to  the  TGT  boimdary  need  not  be  updated  since  all  the  nodes  on  the  layer  just 
inside  the  TGT  boundary  (the  input  ports)  are  zero  for  t>0  (at  t=0  only  one  node  on  the 
layer  just  inside  the  TGT  boimdary  has  the  value  of  1).  The  DBIR  calculations  thus 
require  updates  only  for  the  nodes  between  the  TGT  boundary  and  the  large  grid  boundary 
which  can  result  in  significant  conqiutational  time  savings  if  N»N(^.  The  update 
equations  used  for  the  region  between  the  "inner"  (TGT)  and  the  outer  boundaries  are  the 
firee  space  equations,  since  this  region  is  assumed  to  be  fi:ee  space.  (The  medium  between 
the  two  boundaries  can  be  any  homogenous  medium  extending  to  infinity,  but  the  most 
practical  case  is  the  fi'ee  space.) 

C.  3-D  TGT  RESULTS 

The  overall  procedure  for  3-D  TGT  testing  is  very  similar  to  that  used  for  2-D 
TGT  testhig.  The  mam  difference  is  in  the  amount  of  computer  memory  needed  for  3-D 
TGT  and  "infinite"  grid  moplementation.  Since  3-D  FD-TD  calculations  require  much 
more  computer  memory  because  the  number  of  nodes  in  3-D  is  proportional  to  as 
opposed  to  for  2-D  (N  is  the  number  of  cells  on  each  side  of  the  computationai 
domain)  this  limits  the  grid  size  and  the  duration  of  m^ulse  responses.  We  have  thus 
selected  a  domam  with  only  N=ll  electric  cells  on  each  ade  (electric  cell  faces  are  ahgned 
with  the  domain  boundaries).  Note  that  even  such  a  modest  number  of  cells  per  each  side 
results  in  N^=1331  electric  cells  for  the  domain  vohnne.  Since  each  cell  has  three  field 
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coiDponents,  the  total  niimber  of  electric  field  components  is  3N^=3993.  The  number  of 
magnetic  field  components  is  about  3(N-1)^=3000.  The  total  number  of  unloiowns  is  thus 
about  7000.  3-D  TGT  is  mq)lemented  by  updating  the  electric  field  conq)onents  on  the 
TGT  boimdaiy  based  on  the  time  histories  of  electric  field  conq)onents  at  nodes  just  inside 
the  TGT  bormdary.  The  source  is  a  delta  in5)ulse  of  electric  field  E  field  at  the  center  of 
the  grid,  with  unit-an^Htude  corcponaits  E^,  Ey  and  E^.  The  source  creates  an  outgoing 
spherical  wave.  If  the  grid  and  the  TGT  were  perfect  the  power  within  the  TGT  will  be 
constant  imtil  the  outgoing  wave  has  reached  the  nodes  in  the  centers  of  the  six  TGT 
boundary  planes  (because  these  nodes  are  the  closest  to  the  source  located  at  the  grid 
center.  From  then  on  the  power  within  the  TGT  should  decrease  until  the  outgoing  wave 
has  reached  the  comers  of  the  domain  (the  nodes  fiirthest  away  firom  the  source  at  the 
domain  center).  After  this  time  the  power  within  the  grid,  for  a  perfect  grid  and  TGT, 
should  be  zero.  However,  since  the  grid  is  nor  perfect,  there  will  be  some  residual  power 
within  the  TGT  boundary  even  for  an  infinite  grid.  This  is  the  result  of  the  in^ulse 
"spreading"  as  it  propagates  through  the  FD-TD  grid.  Furthermore,  the  residual  power 
will  include  the  power  "reflected"  off  an  imperfect  TGT.  Again,  just  like  in  2-D,  the 
residual  power  can  be  calculated  as  a  fimction  of  time  for  the  two  cases  infinite  grid  and 
TGT.  The  relative  difference  of  the  two  expressed  in  dB  can  be  thought  of  as  a  measure 
of  performance  of  the  TGT.  This  quantity  can  be  plotted  vs.  time  (with  the  time  ori^ 
defined  as  the  time  when  the  outgoing  wave  first  reaches  the  TGT)  for  various  inpulse 
reponse  durations.  Unfortunately  the  results  for  3-D  TGT  were  not  satisfactory.  The 
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TGT  tests  showed  "late-time  instability",  that  is  the  power  withia  tbe  grid  would  at  first 
decrease  (as  it  should)  but  later  will  start  to  ina-ease  without  bounds,  v^bich  is  incorrect. 
A  saitq>le  result  is  shown  below  (what  is  the  duration  of  the  impulse  response  for  this?) 


Figure  4.14.  Residual  Power  Within  TGT  Showing  Late-Iime  Instability. 

The  non-physical  oscillatory  behaviour  of  the  solution  is  also  evident  firom  the  plot 
of  the  electric  field  at  a  node  inside  the  TGT.  The  field  at  first  converges  to  zero  as  it 
should,  but  after  about  200  steps  it  starts  to  diverge  in  an  oscillatory  manner. 
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Figure  4.15.  Field  Strength  at  a  Node  observation  point  within  boundary. 


For  an  infinite  boundary,  the  residual  power  within  the  TGT  boimdary  decreases  to 
zero  as  it  should.  The  late-time  instability  could  not  be  observed  because  the  available 
PC's  did  not  have  enough  memory  to  accomodate  the  ^e  of  the  required  "infinite"  grid. 
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Figure  4.16.  Residual  Power  within  TGT  for  an  Infinite  Boundary. 
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V.  SUMMARY  AJND  CONCLUSION 


A.  SUMMARY 

In  this  thesis  we  have  tested  the  performance  of  the  TGT  in  1-D,  2-D  and  3-D, 
based  on  the  residual  reflection  off  the  TGT  boundary.  The  power  within  the  TGT 
boundary  was  calculated  at  each  time  step  for  an  "infinite"  grid  (that  is  a  grid  so  large  that 
the  reflection  off  its  boundary  do  not  reach  the  nodes  withiu  the  TGT  boundary  before  the 
calculations  are  completed)  and  for  a  much  smaller  grid  with  the  TGT.  For  1-D  there  were 
no  reflections  off  the  TGT  boundary,  that  is,  the  outgoing  wave  was  "absorbed"  perfectly 
by  the  TGT.  Note  that  in  1-D  there  is  no  in^ulse  spreading.  Therefore,  in  1-D,  the  grid 
and  the  TGT  can  be  "perfect".  This  is  not  the  case  in  2-D  however.  Neither  the  grid  nor 
the  TGT  can  be  perfect,  that  is  there  is  some  pulse  spreading  and  some  "reflection"  off  the 
TGT.  However,  the  power  reflected  off  the  TGT  is  less  than  1%  of  the  residual  power 
due  to  pulse  spreading.  Therefore,  although  the  2-D  TGT  is  not  perfect,  its  performance 
is  excellait  since  the  reflections  off  the  TGT  are  well  below  the  level  of  the  "distortion" 
introduced  by  the  grid  itself  Furthermore,  the  power  reflected  off  the  TGT  boundary 
depends  on  the  diuration  of  the  truncated  impulse  responses:  truncating  uiq>ulse  re^onses 
at  a  later  time  (increasing  the  in5)ulse  response  duration)  reduces  the  reflections  off  the 
TGT.  Finally,  the  3-D  TGT  inplementation  was  not  successful,  because  of  late-time 
mstability  problems. 

B.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  results  of  this  study  show  that  TGT  concept  works  for  1-D  and  2-D  FD-TD 
electromagnetic  field  problems.  The  3-D  TGT  in:q)lementation  was  the  most  difficult, 
because  of  the  corrputer  memory  requirement.  Unfortunately,  it  was  not  successful  due  to 
late-time  instabilities.  More  research  is  needed  to  determine  and  remove  the  causes  of 
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late-time  instabilities  in  3-D.  A  previously  conopleted  thesis  [Wells]  has  demonstrated  the 
effectiveness  of  TGT  for  2-D  static  problems,  v\Me  this  thesis  has  demonstrated  the  TGT 
effectiveness  for  2-D  time-domain  problems.  Another  possibility  for  further  research  is  to 
im5)lement  the  TGT  in  the  frequency  domain,  where  the  time  domain  convolutions  would 
be  replaced  by  sums  of  products  of  coDoplex  numbers. 
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APPENDIX 


A.  DBIR.M 


%%%%%  Program  :  DBIR.M  %%%%% 
%%%%%  This  Program  Generates  TGT  Coefficients  for  2“D  %%%%% 
%%%%%  TMz  fields  (Ez,  Hx,  Hy)  %%%%% 
%%%%%  NOTE:  All  H-fields  have  been  multiplied  by  ZO  I  %%%%% 


clear  all 
clear,  flops (0); 
tic 


N  =  101; 

Nin  =  21; 

Nout  =  Nin  +  2; 


h_dur  =  40; 
vr  =  l/sqrt (2)  ; 
rowsh  =  0.5* (N  -  Nin); 
colsh  =  0.5* (N  -  Nin) ; 
irt  =  rowsh  +  1; 
irb  =  irt  +  Nin  -  1 ; 
id  =  colsh  +  1; 
icr  =  id  +  Nin  -  1; 
nodes__out  =  4*(Nin+l); 
nodes  in  =  4*(Nin-l); 


%  duration  of  impulse  responses 
%  grid  "speed" 

%  Row  "shift”  for  the  "inside"  matrix 
%  Column  "shift”  for  the  "inside"  matrix 
%  Top  row  for  the  "inside"  matrix 
%  Bottom  row  for  the  "inside"  matrix 
%  Leftmost  column  for  the  "inside"  matrix 
%  Rightmost  column  for  the  "inside"  matrix 
%  Number  of  nodes  on  the  boundary 
%  Number  of  nodes  just  inside  the  boundairy 


%%%%%  INITIALIZE  ARRAYS  %%%%%% 

Ez  =  zeros  (N,N);  Hx  =  zeros  (N-2 , N-1)  ;  Hy  zeros  {N-l,N-2)  ; 

Edbir  =  zeros  {nodes_out ,  h__dur)  ;  Dbir  =  zeros  (node s_out ,  nodes_in*h_dur) 
row  pointer  =  zeros  (1 ,  nodes_in)  ;  col_jpointer  =  zeros  ( 1 ,  nodes^in)  ; 
col^start  =  zeros  (1 ,  node s_in)  ;  col_stop  -  zeros  (1 ,  node s__in)  ; 

for  m  =  l:l:nodes_in  %  LCX)P  OVER  INNER  GRID  NODES 

if  (m  <=  Nin) 

row__pointer (m)  =  irt;  col_pointer (m)  -  id  +  m  -  1;  % 

elseif  (m  >  Nin)  &  (m  <=  2*Nin-l) 

row  pointer  (m)  =  irt  +  m  -  Nin;  col^pointer  (m)  icr;  % 

elseif  (m  >  2*Nin-l)  &  (m  <=  3*Nin-2) 

row_pointer  (m)  =  irb;  col__pointer  (m)  =  icr  -  (m- (2*Nin-l)  )  ;  % 

elseif  (m  >  3*Nin-2) 

row_j>ointer  (m)  =  irb  -  (m- {3*Nin-2)  )  ;  col_pointer  (m)  =  id;  % 

end 

col_start(m)  =  (m-l) *h_dur+l ;  col_stop{m)  =  col^start (m)  +  h_dur  - 

end 

pack 


%%%%%%%%%%  DISCRETE  BOUNDARY  IMPULSE  RESPONSE  LOOP  i  %%%%%%%%%%%%% 


%  LOOP  OVER  INNER  GRID  NODES  as  sources 
for  m  =  l:l:nodes_in 

Ez  zeros  (N,N);  Hx  =  zeros  (N-2  ,N-1)  ;  Hy  =  zeros  (N-l,N-2); 


TOP  SIDE 
RIGHT  SIDE 
BOTTOM  SIDE 
LEFT  SIDE 
1; 
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Ez  {row_pointer  (m)  ,  col_jpointer  (m)  )  =  1  ; 


%  Delta  pulse  excitation  along  the.  top 
half  of  the  boundary 

%%%%%%  TIME  LOOP  %%%%%% 
for  k  =  l:l:h_dur 

Etop  =  (Ez (irt-l, icl>l : icr+1) ) ' ; 

Bright  =  Ez (irt :irb+l , icr+1) ; 

Ebot  =  (fliplr (Ez (irb+1, icl“l :icr) ) ) ' ; 

Eleft  =  flipud (Ez (irt : irb, icl“l) ) ; 

Edbir(:,k)  =  [Etop;  Bright;  Ebot;  Eleft]; 

Hx  =  Hx  +  vr*(  Ez(2:N-l,l:N-l)  -  Ez (2 :N-1, 2 :N) ) ; 

Hy  =  Hy  +  vr* ( -Ez (1 :N-1 , 2 :N-1 )  +  Ez (2 :N, 2 :N-1) ) ; 

Ez  (2  :N-1, 2  :N-1)  :=Ez  (2  :N-1, 2  :N-1)  +vr*  (-Hy{l  :N~2,  : )  +Hy  (2  :N-1,  : ) 

+Hx( : ,l:N-2) -Hx( :,2:N-1)); 

Ez  (irt :  irb,  id  :  icr)  =  zeros  (Nin,  Nin)  ;  %  Reset  Ez  inside 
end 

%%%%%%  EK:j  OF  TIME  LOOP  %%%%%% 

Dbir(;,  col_start  (tn)  :  col_stop  (m)  )  =  Edbir; 

end 

%%%%%%  END  OF  SOURCE  LOOP  %%%%%% 

MFlops  =  flops/1000000 

toe 

save  dbir_40Tn  Dbir 

save  par_4  0in  Nin  h_dur  vr 


B.  TGT  TM.M 


%%%%%%  Program  :  TGT_TM.M  %%%%%% 
%%%%%%  FD-TD  for  2D,  TMz  fields  (Ez,  Hx,  Hy)  ,  cartesian  cordinates  %%%%%% 
%%%%%%  NOTE:  H-fields  have  been  multiplied  by  ZO  !  %%%%%% 
%%%%%%  Dirichlet  b.c.  or  OPEN  boundary  simulated  via  DBIR  %%%%%% 


load  dbir_40m 
load  par_40m 

%%%%%  INPUT  PARAMETERS  %%%%% 

N  =  Nin  +  2; 

nodes  out  =  4*(Nin+l);  %  Number  of  nodes  on  the  boundary 

nodes_in  =  4*(Nin-l);  %  Number  of  nodes  just  inside  the  boundary 

%  Time  Samples 
N_samples  =  200; 

%%%%%  ARRAYS  INITIALIZED  %%%%% 

Ez  =  zeros (N) ;  Hx  =  zeros (N- 2 ,N-1) ;  Hy  =  zeros (N-l,N-2) ; 

Ebound  =  zeros  (nodes_in*h_dur,  1)  ;  E_dbir  =  zeros  (nodes_out,  1)  ; 
time_pointer  ==  zeros  {nodes_in,  1)  ; 

%%%%%  BOUNDARY  POINTERS  FOR  THE  TIME  HISTORIES  OF  NODES  JUST  INSIDE  %%%%% 

for  i  =  l:nodes_in,  time^pointer (i)  *  (i-l)*h_dur  +  1;  end 

%%%%%  BOUNDARY  POINTERS  FOR  SPIRAL  TO  RASTER  LABELING  CCMIVERSION  %%%%% 
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%%%%%  (can  be  done  only  once  and  storedlli)  %%%%% 

for  n  =  1 : 1  :nodes__out 
if  n  <=  N 

rowjpointer  (n)  =  1;  col_jpointer  (n)  =  n; 

elseif  (n  >  N)  &  (n  <=  2*N-1) 

row_j>ointer  (n)  =  n  -  N  +  1;  col_pointer  (n)  =  N; 
elseif  (n  >  2*N-1)  &  (n  <=  3*N-2) 

rowjpointer  (n)  =  N;  col_pointer  (n)  =  3*N  -  1  -  n; 

elseif  (n  >  3*N-2) 

roWjPointer  (n)  ~  4*N  -  2  -  n;  col_pointer (n)  =  1; 
end 

end 


%%%%  Computation  of  Residual  Power  for  grid  %%%% 

%%%%  since  wavefront  touching  the  tgt  boundary  %%%% 

if  (k>=ceil(N/2)  &  k<=ceil (N/2) +140) 

m:=m+l ; 
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6  (m)  =0  ; 
for  a=l:N 
for  b=:l:N 

& (m) -s (m) +Ez (a, b) ^2 ; 
end 

end 


end 


end 

save  s_tgt  s 
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